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CHAPTER  I 

INTRODUCTION 

The  transonic  flow  about  an  ax i symmetric  jet  engine 
nacelle  can  create  pressure  forces  which  may  contribute  sig- 
nificantly to  the  total  drag  of  the  aircraft.  In  wind  tunnel 
testing  of  such  configurations,  it  is  desirable  to  have 
analytical  methods  to  predict  the  pressure  distribution  over 
the  nacelle.  With  reliable  prediction  methods,  areas  of 
specific  interest  in  the  Mach  number /Reynolds  number  test 
envelope  may  be  investigated.  Such  methods  are  also  of  value 
for  the  placement  and  sizing  of  wind  tunnel  instrumentation 
and  an  evaluation  of  wall  interference  effects. 

Particular  interest  has  centered  around  the  aft 

portions  of  the  nacelle  under  the  general  topic  of  nozzle- 

afterbody  (NAB)  flow  problems.  The  combination  of  internal 

and  external  flows,  with  their  viscous  properties,  and  the 

♦ 

possibility  of  a shock  wave  on  the  afterbody  creates  a flow 
that  is  complex  and  highly  interactive.  The  jet  plume  it- 
self, considered  in  its  relationship  with  the  external  flow, 
has  the  dual  characteristics  of  a solid  body,  causing  turning 
of  the  external  flow,  and  of  a shear  layer,  which  accelerates 
the  external  flow  due  to  the  higher  jet  velocity,  in  general. 
Any  prediction  must  include,  in  some  way,  all  of  the  above- 
mentioned  characteristics  if  it  is  to  be  accurate. 
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In  order  to  render  this  complex  flow  manageable  to 
analysis,  an  approximate  model,  as  shown  in  Figure  1,  is 
assumed.  Most  methods  that  have  addressed  themselves  to  the 
solution  of  the  NAB  problem  use  a similar  type  of  flow 
model  (1)^,  differing  primarily  in  the  way  that  the  mixing 
region  is  handled.  Using  this  model,  the  usual  approach  is 
to  use  a viscous/inviscid  iterative  technique,  in  which  the 
body  shape  is  modified  by  the  boundary-layer  displacement 
thickness  calculated  along  the  length  of  the  body.  An  in- 
viscid  calculation  of  the  external  flow  is  then  carried  out 
on  the  new  body  shape,  and  the  viscous/inviscid  iteration  is 
continued  until  convergence.  The  plume  shape  may  be  modified 
in  some  way  due  to  the  effects  of  flow  entrainment  and  changes 
in  the  pressures  of  the  local  inviscid  flow.  If  the  nacelle 
IS  mounted  on  a solid  sting,  then  no  entrainment  analysis  is 
required  and  the  viscous/inviscid  method  is  directly  appli- 
cable to  the  nacelle/sting  combination. 

Numerous  recent  efforts  have  used  a viscous- 
inviscid  iterative  method  to  predict  pressure  distributions 
with  various  degrees  of  success.  Grossman  and  Melnik  (2) 
were  among  the  first  to  publish  results  of  a combined  viscous 
inviscid  iterative  computer  program,  but  they  treated  the  net 
plume  as  a solid  body.  Further  efforts  at  Grumman  (3)  took 
the  entrainment  effect  of  the  plume  in-to  account  through  the 


^Numbers  in  parentheses  refer  to  similarly  numbered 
references  in  the  Bibliography. 
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Figure  1 . Schematic  of  the  various  flow  regions  involved  in  the 
modeling  of  the  afterbody  problem. 
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use  of  a correction  to  the  plume  shape  based  on  Green's  inte- 
gral method  (4).  This  method  is  still  the  basis  of  the  plume 
entrainment  model,  according  to  the  most  recent  publica- 
tion (5) , but  current  efforts  have  shifted  to  the  case  of 
afterbodies  with  flow  separation  occurring  on  the  afterbody, 
with  reattachment  located  on  the  solid  supporting  sting  or 
plume.  Presz's  method  (6)  for  predicting  and  calculating  the 
separated  region  is  used.  It  is  significant  to  note  that 
even  in  cases  where  the  flow  is  substantially  attached, 
Grossman  and  Melnik  felt  that  the  Green's  integral  method 
"...  [did]  not  take  into  account  the  complete  mixing  of  the 
high  velocity  jet  with  the  external  stream.  Nonetheless, 
this  procedure  [did] , at  least  qualitatively,  produce  a 
necking-down  of  the  wake  displacement  surface  and  [indicated] 
some  jet-entrainment  effects."  (3) 

Keith,  et  al.  (7)  addressed  the  more  complicated 
problem  of  the  transonic  flow  about  a two-dimensional  or  axi- 
symmetric  high-bypass-ratio  fan  duct  nacelle  with  appropriate 
inlet  and  exhaust  flows.  Their  solution  technique  was  a 
Streamtube  Curvature  Analysis,  in  which  the  streamline 
coordinates  were  updated  based  on  new  coordinates  required  by 
the  continuity  equation  and  the  effect  of  the  resulting 
change  in  streamline  curvature  as  computed  from  the  integrated 
momentum  equation.  This  inviscid  solution  technique  was  com- 
bined with  a boundary-layer  analysis,  which,  however,  was  used 
only  over  the  solid  body  to  provide  a new  effective  body  for 
the  inviscid  analysis.  This  work  also  included  calculations 
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for  incipient  separation,  but  no  detailed  analysis  of  the 
flow  field  from  separation  to  recompression  was  carried  out. 

In  addition,  any  changes  in  jet  plume  shape  itself  were 
neglected. 

A "wake-body"  configuration  for  flow  over  an  after- 
body with  a blunt  base  and  no  plume  was  used  by  Rom  and 
Bober  (8)  in  their  iterative  solution  technique.  At  the  end 
of  the  boattail  a conical  body  was  assumed  over  which  the 
flow  was  calculated  in  the  regular  viscous/inviscid  iterative 
manner  until  the  boundary  layer  solution  began  to  diverge. 

The  minimum  radius  was  then  extended  as  a solid-body  cylinder 
downstream.  For  uniqueness  of  solution,  however,  an  experi- 
mental pressure  at  the  separation  point  was  required,  as  the 
shape  of  the  "wake-body"  was  varied  until  this  pressure  was 
matched.  The  "wake-body"  concept  will  be  used  in  this  paper 
for  one  configuration  similar  to  those  configurations  studied 
by  Rom  and  Bober.  In  a broader  sense,  however,  it  can  be 
said  that  the  wake-body  approach  is  basic  to  the  modeling  of 
plume  or  wake  flow,  whether  the  flow  is  separated  or  at- 
tached, when  the  viscous/inviscid  iterative  method  of  solu- 
tion is  utilized. 

It  is  striking  to  see  the  advances  made  in  recent 
years , both  from  the  recognition  of  the  similarity  of  ap- 
proach, vis-a-vis  the  iterative  technique,  to  the  application 
of  the  particular  differences  in  the  detailed  problems  con- 
sidered by  each  effort.  Chow,  Bober,  and  Anderson  (9)  pointed 
out  the  inadequacy  of  small  disturbance  theory  for  use  in  the 
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inviscid  portion  of  the  iteration.  Their  predictions  of 
attached  flows  over  a boattail  with  a cylindrical  sting  were 
improved  when  the  full  transonic  potential  equations  were 
used.  They  were  also  among  the  first  to  acknowledge  the 
necessity  of  "relaxing"  the  boundary  layer  correction  as  a 
means  of  attaining  convergence  of  the  iterative  technique. 
Bower  (10)  applied  the  iterative  technique  to  the  problem  of 
ax i symmetric  or  two-dimensional  diffusers.  His  method  is 
interesting  in  that  he  applied  Strong  Interaction  Theory  once 
his  solution  indicated  separation,  that  is  to  say,  when  his 
viscous  solution  began  to  diverge,  he  would  begin  calculating 
the  flow  step-by-step  downstream  using  a one-dimensional  de- 
scription of  the  inviscid  region  solved  simultaneously  with 
the  integral  equations  for  the  turbulent  compressible  viscous 
layer.  Calarese  (11)  approached  the  jet  plume  displacement 
phenomenon  by  means  of  varying  the  location  of  a cone  frustrum 
on  the  NAB  support  sting  according  to  a correlation  given  by 
McDonald  and  Hughes  (12) . Any  separated  regions  were  calcu- 
lated using  a Korst  base  flow  analysis.  Calarese 's  method, 
however,  was  limited  to  subsonic  free-stream  flows. 

In  the  work  of  Moulden,  Wu , and  Spring  (13)  and 
the  more  comprehensive  work  of  Wu , Moulden,  and  Uchiyaraa  (14) , 
excellent  descriptions  of  the  interaction  of  the  flow  com- 
ponents about  a missile-type  configuration  are  given.  The 
approach  to  the  solution  of  the  entire  flow  field  about  a 
missile  with  and  without  fins  is  the  viscous/inviscid  itera- 
tion, but  the  inviscid  portion  of  the  analysis  is  somewhat 
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limited  in  that  the  transonic  small  perturbation  equation  has 
been  used  throughout.  The  attached  boundary  layer  flow  is 
handled  in  a straightforward  manner  employing  the  method  of 
Nash  (15)  , while  the  inviscid  plume  shape  is  calculated  with 
a step-by-step  method  of  characteristics  solution  taking  into 
account  the  static  pressure  variation  due  to  the  interaction. 
Concerning  itself  with  missile  configurations,  however,  the 
complete  method  assumes  a strong  interaction  between  the  jet 
plume  and  the  external  flow,  and  the  analysis  assumes  a size- 
able separated  area  from  the  afterbody  to  a point  of  reattach- 
ment, or  "confluence,”  on  the  plume.  This  region  is  treated 
with  a Korst  component  analysis. 

Chronologically,  the  past  year  has  witnessed  the 
culmination  of  what  may  be  regarded  as  the  first  stage  of  the 
viscous/inviscid  iterative  solution  of  transonic  nozzle  after- 
body flows,  including  the  effects  of  shock  boundary  layer 
interaction,  plume  and  external  flow  interaction,  and  after- 
body flow  separation  and  reattachment.  This  was  demonstrated 
at  the  Propulsion  interactions  Workshop  held  at  Langley 
Research  Center  in  May  1976  (16) , where  no  less  than  five 
papers  were  presented  dealing  with  the  viscous/inviscid  iter- 
ative method,  and  a number  of  other  papers  were  delivered 
treating  the  analyses  of  individual  components  of  the  flows. 
The  most  recent  papers  also  indicate  a continued  and  wide- 
spread interest  and  confidence  in  iterative  techniques. 

Reubush  and  Putnam  (17)  and  the  continuing  efforts  related  to 
this  work,  Putnam  and  Abeyounis  (18)  and  Wilmoth  (19) , have 
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applied  the  iterative  method  to  flows  with  separation,  accord- 
ing to  the  Presz  method,  but  have  not  included  any  plume  en- 
trainment model  in  their  analysis.  Cosner  and  Bower  (20) 
have  also  concentrated  on  the  separated  case,  using  an 
"empirical  bridging"  technique  to  carry  the  flow  analysis 
past  the  point  of  reattachment  to  a point  where  boundary  layer 
calculations  may  once  again  be  used.  Plume  entrainment  is 
modeled  by  treating  the  inviscid  plume  boundary  as  a non- 
adiabatic  moving  wall.  The  present  author's  most  recent 
utilization  of  the  viscous/inviscid  iteration  (21)  included  a 
crude  model  for  plume  entrainment  based  on  mass  conservation 
considerations.  Finally,  a most  interesting  approach,  and 
one  that  may  be  representative  of  future  efforts,  was  pro- 
vided by  Holst  (22),  who  used  the  iteration  technique  com- 
bined with  solution  of  the  Navier-Stokes  equations  in  the 
afterbody/plume  region.  A summary  of  the  above  information  is 
presented  in  Table  I. 

The  purpose  of  this  thesis  will  be  to  provide  a 
viscous/inviscid  iterative  technique  valid  for  the  solution 
of  transonic,  attached  flows  over  typical  jet  engine  nacelles, 
including  a reliable  and  efficient  method  of  accounting  for 
the  effect  of  jet  plume  entrainment  on  the  afterbody  pressure 
distribution.  Accordingly,  Chapter  II  will  describe  the 
development  of  the  combined  viscous  and  inviscid  digital  com- 
puter program  and  Chapter  III  will  derive  the  equations,  based 
on  a Chapman-Korst  mixing  analysis,  necessary  for  the  imple- 
mentation of  the  plume  entrainment  model.  Results  of  the 
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SUMMARY  OF  PREVIOUS  VISCOUS/TNVISCID  ITERATIVE  EFFORTS 


Name  and  Reference 

Plume 

Entrainment 

Separated 

Flows 

Comments 

Grossman  & Melnik  (2) 

No 

No 

Treated  plumes  as  solid  body 

Grossman  & Melnik  (3) 

Yes 

No 

Green's  Integral  Method  for  mixing  ' 

Yaeger  (5) 

Yes 

Yes 

Presz's  Method  for  separated  flows 

Keith,  et  al . (7) 

No 

No 

Streamtube  curvature  analyses  of 
( bypass  jet  engine 

Rom  and  Bober  (8) 

No 

Yes 

Wake  body  for  separated  region 

Chow,  Bober,  and 

Anderson  (9) 

No 

No 

"Relaxed"  boundary  layer 

Bower  (10) 

Cosner  & Bower  (20) 

NO 

Yes 

Strong  interaction  theory  for 
separated  regions 

Calarese  (11) 

No 

Yes 

One-stream  Chapman-Korst  analysis 
for  separated  region 

Wu , Moulden , and 

Uchiyama  (14> 

Yes 

Yes 

Separated  flow  assumed 

Reubush  & Putnam  (17) 
Putnam  & Abeyounis (18) 

Yes 

No 

Presz's  Method  for  separated  flows 

Yaros  (21) 

Yes 

No 

Plume  entrainment  effect  based  on 
mass  conservation 

Holst  (22) 

No 

Yes 

Navier-Stokes  equation  used  in  after- 
body and  plume  region 
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combined  program,  compared  with  data  from  typical  practical 
geometric  configurations,  will  be  presented  in  Chapter  IV, 
with  Chapter  V presenting  conclusions  and  recommendations 
based  on  these  efforts. 
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CHAPTER  II 

THE  VISCOUS/INVISCID  ITERATIVE  TECHNIQUE 
AND  THE  COMBINED  PROGRAM 

Introductory  Remarks 

This  section  will  describe  the  physical  model  of 
the  transonic  viscous  flow  over  a typical  nozzle/afterbody 
combination,  emphasizing  the  case  where  the  flow  is  largely 
attached.  Further  simplifications  of  the  model  to  enhance 
the  mathematical  resolution  of  the  complex  flow  field  will 
then  be  presented,  followed  by  an  explanation  of  the  viscous/ 
inviscid  iterative  technique  as  applied  in  this  thesis.  A 
summary  of  the  component  inviscid,  viscous,  and  pltime  digital 
computer  programs  will  indicate  the  advantages  and  limita- 
tions of  this  method,  and  the  iterative  technique  as  carried 
out  operationally  with  the  combined  program  will  be  described. 

Viscous/Inviscid  Iterative  Approach 

As  indicated  previously,  the  physical  model  applied 
to  transonic  flows  over  an  afterbody  and  plume  combination 
usually  takes  the  form  shown  in  Figure  1,  page  3.  The  flow 
over  the  afterbody  and,  for  that  matter,  the  flow  over  the 
body  preceding  the  afterbody  are  broken  down  into  their  in- 
viscid and  viscous  components,  a result  of  the  application  of 
classical  boundary  layer  principles.  This  simplification  of 
the  flow  over  the  body  is  justifiable  if  no  local  flow 
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separations,  due  to  strong  shock-boundary  layer  interactions 
and/or  excessive  changes  in  body  shape,  are  present.  The 
flow  within  the  nozzle  is  approximated  in  an  identical  manner, 
with  similar  restrictions  on  local  flow  separation-  The 
resulting  flow  at  the  end  of  the  nozzle/afterbody  is  still 
within  the  realm  of  the  boundary  layer  concept,  that  is, 
there  exists  two  boundary  layers  external  to  which  are  their 
respective,  essentially  inviscid,  external  and  internal  flows. 

With  no  viscous  effects,  for  example  if  the  internal 
and  external  flows  were  fully  inviscid  and  uniform  at  the  end 
of  the  nozzle/afterbody,  the  subsequent  flow  development  would 
follow  the  line  of  the  inviscid  jet  plume,  denoted  in  Figure 
1,  page  3,  as  the  inviscid  jet  plinne  boundary.  This  boundary 
is  a function  of  the  ratio  of  the  pressure  of  the  jet  and  the 
pressure  of  the  external  stream  as  well  as  the  thermodynamic 
properties  of  both  streams.  Under  this  inviscid  assumption, 
the  jet  plume  boundary  would  act  as  a solid  obstruction  to 
the  external  flow,  causing  it  to  turn  more  rapidly  in  the 
neighborhood  of  the  nozzle/af terbody  endpoint,  thus  increasing 
the  pressure  of  the  afterbody-  This  phenomenon  is  characterized 
as  the  "displacement  effect"  of  the  plume  shape. 

In  the  actual,  physical  case  there  exists  a strong 
viscous  interaction  between  the  internal  and  external  streams 
after  the  nozzle/afterbody  endpoint.  There  is  substantial 
turbulent  mixing  between  the  two  streams  with  mass  and  momen- 
tum transfer  across  the  inviscid  jet  boundary.  This 
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phenomenon,  known  as  the  "entrainment  effect"  of  the  jet 
plume,  will  change  the  effective  shape  of  the  plume  as  seen 
by  the  external  inviscid  flow. 

It  is  for  these  reasons  that  the  flow  model  shown 
in  Figure  2 is  assumed  for  purposes  of  applying  the  viscous/ 
inviscid  iterative  technique.  The  technique  is  restricted  to 
afterbodies  with  sharp  trailing  edges,  i.e.  no  substantial 
base  areas  are  permitted.  Basically,  this  model  reflects  the 
assumption  that  the  plume  boundary  will  be  treated  as  a solid 
boundary  in  the  first  iterative  computation.  As  a qualitative 
example,  consider  the  afterbody  and  original  plume  boundary 
as  shown  in  Figure  2.  Application  of  the  pressure  distribu- 
tion over  this  body /plume  combination  to  the  viscous  boundary 
layer  solution  would  result  in  a classic  viscous/inviscid 
gQl'^'tion  as  shown  by  the  solid  lines  denoting  the  displacement 
thickness  of  the  external  flow.  If  the  viscous  mixing  between 
the  plume  and  external  flow  were  included,  however,  the 
original  plume  shape  would  be  modified  as  shown  by  the  dashed 
line,  as  would  be  the  resulting  effective  body  as  determined 
by  the  boundary  layer  displacement  thickness  and  modified 
plume  shape . 

One  is  then  led  naturally  to  the  concept  of  the 
viscous/inviscid  iterative  technique.  With  the  new  effective 
body  and  effective  plume  shape  the  inviscid  calculation  is 
repeated.  If  the  pressure  distribution  over  the  body  and 
plume  agrees  with  the  previous  pressure  distribution,  the 
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Figure  2.  Schematic  of  the  modified  flow  regions  involved  in  the 
modeling  of  the  afterbody  problem. 
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solution  is  considered  converged.  If  not,  the  viscous/in- 
viscid  iteration  is  continued. 

One  advantage  to  this  approach  is  that  the  viscous/ 
inviscid  iterative  operation  may  be  developed  for  flows  over 
what  are  effectively  solid  bodies,  and  the  adjustment  to  the 
effective  plume  shape,  which  will  be  developed  in  Section  III, 
may  be  included  when  necessary.  Thus,  flows  over  afterbodies 
with  solid  support  stings  may  be  treated  directly  by  this 
method  by  simply  not  applying  the  plume  entrainment  adjust- 
ments . 

Inviscid  Digital  Computer  Program 

The  use  of  a reliable  and  accurate  inviscid  digital 
computer  program  in  the  viscous/inviscid  iterative  method  is 
of  prime  importance,  in  that  it  is  only  through  this  segment 
of  the  iteration  that  any  upstream  influences  of  perturbations 
of  the  flow  can  be  accounted  for,  as  the  viscous  segment  is 
essentially  parabolic  in  the  longitudinal  direction.  Certain 
requirements  peculiar  to  the  nature  of  the  transonic  after- 
body flow  problem  demand  an  inviscid  digital  computer  program 
with  the  following  characteristics: 

1.  The  program  must  be  fully  transonic. 

2.  The  program  must  be  axisymmetric . 

3.  Small  perturbation  approximations  are  valid 
for  restricted  body  shapes  and  may  not  be 
satisfactory  for  the  afterbody  geometries 
considered. 
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4 . Arbitrary  effective  body  shapes  must  be 
acceptable  to  the  program, 

5.  Computation  time  must  not  be  excessive. 

In  early  applications  of  the  iterative  method  (23) , 
three  digital  computer  programs  were  used  and  compared,  two 
of  which  proved  to  be  unsatisfactory.  Of  these  two,  the 
first,  the  Douglas-Neumann  Potential  Flow  Program,  displayed 
favorable  accuracy  and  reliability,  but  it  was  limited  to 
wholly  subsonic  flow.  The  second  program  was  the  Pratt  and 
Whitney  Transonic  Program,  which  was  based  on  the  small  per- 
turbation form  of  the  transonic  flow  equations.  In  addition 
to  this  restriction,  it  displayed  other  shortcomings  related 
to  extreme  sensitivity  to  the  effective  body  shape  derivatives 
and  a chronic  inability  to  provide  reliable  solutions  for 
supersonic  portions  of  the  flow.  By  far  the  most  acceptable 
inviscid  computer  program  examined  was  the  third,  the  South- 
Jameson  Relaxation  Solution  for  Inviscid  Axisymmetric  Flow, 
RAXBOD  (24).  RAXBOD  is  a finite-difference  relaxation  solu- 
tion of  the  full  transonic  potential  equation  with  exact 
boundary  conditions.  The  program  is  applicable  only  to  axi- 
symmetric bodies  in  steady  flow.  These  characteristics  of 
RAXBOD  satisfy  the  first  three  requirements  stated  previously. 

In  addition,  RAXBOD  contains  some  other  special 
features  which  enhance  its  use  for  flow  calculations  about 
afterbody  and  sting  or  plume  combinations,  as  well  as  for 
flow  about  the  entire  axisymmetric  body.  Originally  conceived 
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as  a calculation  method  for  blunt  axisymmetric  bodies,  RAXBOD 
uses  two  different  coordinate  systems  along  the  length  of  the 
body,  as  shown  in  Figure  3.  Based  on  the  experience  of 
Lipnitskii  and  Liftshits  (25),  South  and  Jameson  chose  to  use 
a conformal  body  normal  coordinate  system  in  the  body  nose 
region.  This  coordinate  system  supplies  the  fine  and  well 
balanced  grid  distribution  necessary  in  this  region  of  highly 
changeable  flow.  Surprisingly,  this  coordinate  system  also 
works  well  for  pointed-nose  bodies  by  assuming  the  body  angle 
to  be  90  degrees  at  the  nose  and  supplying  the  other  normal 
coordinate  axes  at  their  true  body-normal  angles.  According 
to  South  and  Jameson,  "This  procedure  leaves  an  embarrassing 
gap  in  the  coordinates  in  the  region  ahead  of  the  nose,  but 
the  calculation  proceeded  without  difficulty."  (24)  The 
present  author  has  also  experienced  little  difficulty  with 
calculations  in  the  region  of  a pointed  nose,  typical  of  an 
equivalent  body  of  revolution  rather  than  a nacelle,  provid- 
ing that  the  shape  is  convex,  in  spite  of  the  fact  that  the 
difference  in  orientation  of  the  first  two  normal  axes  is 
often  75  degrees  or  more. 

At  some  distance  down  the  body,  usually,  but  not 
always,  the  point  where  the  body  obtains  a horizontal  slope, 
the  sheared  cylindrical  system  of  coordinates  is  introduced . 
In  this  system,  the  coordinate  axes  extending  out  into  the 
flow  remain  perpendicular  to  the  axis  of  symmetry  of  the 
body,  while  the  tangential  coordinate  axes  remain  their  re- 
spective distances  from  the  body  at  each  grid  intersection. 
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Figure  3*  Two  coordinate  systems  used  by  RAXBOD. 
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thus  attaining  a "sheared"  configuration  matching  the  body 
shape.  This  coordinate  system  was  expressly  used  to  avoid 
the  problems  a body-normal  system  would  experience  in  areas 
of  body  concavity,  where  two  or  more  body-normal  coordinates 
would  tend  to  intersect  at  large  distances  from  the  body, 
creating  computational  difficulties.  Thus,  with  the  excep- 
tion of  some  points  in  the  nose  region  where  the  body-normal 
system  is  used,  RAXBOD  is  capable  of  flow  calculation  over 
many  varied  body  shapes  and  is  thus  satisfactory  from  the 
standpoint  of  the  fourth  requirement. 

Another  feature  of  RAXBOD  which  lends  itself  to 
afterbody  configurations  with  highly  sloping  shapes  is  the 
"Rotated  Difference  Scheme,"  which  is  applied  at  supersonic 
points  when  both  velocity  components,  normal  and  tangential 
to  the  local  coordinates,  are  subsonic.  Although  South  and 
Jameson  point  out  that  the  Rotated  Difference  Scheme  is  not 
an  exact  transformation  to  an  intrinsic  coordinate  system,  it 
does  provide  a redefinition  of  the  computational  "molecule" 
to  include  new  points  when  the  flow  is  not  well  aligned  with 
the  computational  coordinates. 

The  RAXBOD  Progreim  thus  solves  the  equation 


1 ^ 
H 3? 


1 

H 35 

L 


T UV  3^4. 


+ 


19 


AEDC-TH-77-106 


in  the  orthogonal  curvilinear  coordinate  system,  where 


and 


U = cos  0 + i II 


V — “ sin  0 + -jI  ■ 
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In  the  sheared  cylindrical  system,  the  equation 
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The  boundary  conditions  are  applied  at  infinity 
in  the  forms 


(j>  0 as  n 

in  both  coordinate  systems,  and,  additionally. 
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e -►  0 as  ^ “ 

when  the  sheared  cylindrical  coordinates  are  used.  The 
boundary  conditions  at  the  body  are 

V = 0 -S'  = sin  9 


for  the  orthogonal  curvilinear  coordinates,  and 


S'  - “bod  - 0 * = 


■^bod 
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for  the  sheared  cylindrical  system.  The  computational  pro- 
cedure of  South  and  Jameson  also  allows  for  coordinate 
stretching  in  the  body  normal  direction  near  the  body,  and  an 
independent  stretching  in  the  tangential  direction  in  regions 
where  large  gradients  of  the  perturbation  potential  are  ex- 
pected to  occur,  e.g.  the  nose  aind  afterbody  regions. 

The  Rotated  Difference  Scheme  is  applied  in  areas 
of  supersonic  flow,  i.e. 


The  basis  behind  the  scheme  is  to  use  upwind  differences  for 
the  second  derivatives  contributing  to  cjigg  and  central 
differences  for  the  contributions  to  'tjgjji  where  these 
derivatives  are  in  the  directions  tangential  to  and  normal  to 
the  velocity  vector,  respectively.  The  derivatives  for  the 
orthogonal  curvilinear  system  become 
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and 
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The  corresponding  equations  for  the  sheared  cylin- 
drical system  are 
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The  upwind  differences  for  the  (fgg  derivatives  are  taken  in 
the  quadrant  in  which  the  velocity  vector  is  located. 

Computation  time,  requirement  5,  is  also  a factor  in 
any  relaxation  method.  South  and  Jameson  used  the  steady 
transonic  equations  with  spatial  relaxation  as  opposed  to  the 
time-asymptotic  method,  based  on  findings  by  Yoshihara  (26) 
that  the  latter  method  was  more  costly  in  computer  time  and 
storage.  To  further  reduce  computation  time,  RAXBOD  begins 
calculations  with  a coarse  computational  grid,  for  example 
25  X 25  points,  and  upon  convergence  in  this  grid  system, 
halves  the  mesh  size  and  again  calculates  to  convergence. 
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This  "mesh  halving”  may  be  repeated,  bringing  the  final  mesh 
to  a 97  K 97  point  configuration,  at  a considerable  saving, 
according  to  South  and  Jameson,  of  computation  time  if  the 
97  X 97  grid  were  used  throughout  (Table  II) . The  present 
author's  experience,  based  on  the  two  mesh  halvings  described 
above,  is  shown  in  Table  III.  While  RAXBOD  may  not  be  con- 
sidered an  extremely  fast  program,  the  computation  times  in- 
volved are  reasonable  when  the  complexity  and  sensitivity  of 
the  calculation,  as  well  as  the  size  of  the  computational  grid, 
are  considered. 

The  RAXBOD  program  was  originally  designed  to  cal- 
culate flow  over  a solid  body,  the  shape  of  which  was  well 
defined,  and  as  such  required  the  body  shape  to  be  input  as 
an  analytic  function  or  combinations  of  such  functions,  such 
as  circular  arcs,  cylinders,  or  cones.  The  first  and  second 
derivatives  at  each  longitudinal  grid  location  were  then  inte- 
grated along  the  length  of  the  body  to  provide  the  body  shape 
for  computation.  As  the  non-analytic  nature  of  the  effective 
body  shape  during  the  course  of  a viscous/inviscid  iteration 
did  not  lend  itself  to  such  specification,  the  program  was 
modified  to  accept  body  coordinates  as  input.  The  required 
derivatives  calculated  from  neighboring  points  were  provided 
to  the  program  for  subsequent  integration  using  a curve 
fitting  technique  based  on  the  work  of  Akima  (27) . It  is 
worth  noting  that  the  most  recent  version  of  RAXBOD  (28)  has 
been  modified  to  accept  body  coordinates  directly. 
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TABLE  II 

CONVERGENCE  HISTORY  WITH  AND  WITHOUT  MESH 
HALVING  (ADAPTED  FROM  REF.  24) 


Computational  Grid 

Cycles 

Adk  a 

^max 

Time^ 

With  Mesh  Halving 

25  X 25 

36 

4.9  X 10"^ 

6.6 

49  X 49 

21 

4.7  X 10"^ 

15.3 

97  X 97 

53 

9.8  X 10“® 

154.0 

No  Mesh  Halving 

97  X 97 

280 

1.0  X 10“^ 

812.0 

= max 
^itiax  ' 

^1,3  ^1,3 

, n = number 

of  cycles  and 

. . is  the  transonic  disturbance  potential  at 

t z 

grid  point 

(i , j)  . 


i- 

CDC  6600  CPU  time,  seconds. 
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TABLE  III 

TYPICAL  COMPUTATION  TIMES  USING  MESH  HALVING 


Mach 

Number 

(25x25) 

Cycles 

(49x49) 

(97x97) 

A A a 

A d) 

^max 

Time^ 

0.6 

5 

16 

30 

1.0  X 10"^ 

48.8 

0.8 

5 

19 

30 

44.2 

0.9 

6 

26 

31 

47.1 

0.95 

16 

38 

48 

68.7 

1,05 

18 

26 

101 

157.2 

1.1 

20 

25 

121 

184.7 

1.2 

23 

30 

119 

202.0 

1.4 

42 

68 

237 

435.8 

^A(ti  defined  as  in  Table  II. 


W 

“iBM  370  mod  165  time,  seconds. 
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viscous  Digital  Computer  Program 

The  viscous  portion  of  the  viscous/inviscid  iter- 
ative method  is  no  less  important  than  the  inviscid  portion. 
Although  there  exist  literally  dozens  of  "boundary  layer 
programs"  capable  of  calculation  of  the  viscous  segment  of 
the  flow,  the  requirements  of  the  iteration  demand  certain 
characteristics  of  the  program  which  may  be  summed  up  below: 

1.  the  prograun  must  be  axisymmetric , 

2.  the  program  must  be  able  to  calculate 
through  strong  adverse  pressure  gradients 
which  are  generated  by  the  inviscid  program 
in  the  early  iterations,  and 

3.  computation  time  must  not  be  excessive. 

The  five  viscous  digital  computer  programs  examined 
were  essentially  different,  as  a representative  cross-section 
of  solution  methods  was  desired.  For  purposes  of  comparison, 
a hypothetical  afterbody/plume  shape  was  assumed  (Figure  4) , 
and  each  program  was  run  under  the  conditions  of  zero  longi- 
tudinal pressure  gradient.  The  body  shape  was  represented  in 
each  program  by  discrete  points,  and  thus  the  sharp  turns  in 
Figure  4 were  not  really  present.  The  five  viscous  methods 
examined  were  those  of  Kuhn  and  Nielsen  (29)  , Tucker  (30) , 
Bartz  (31),  Mayne  (32),  and  Whitfield  (33).  Of  these,  the 
method  of  Kuhn  and  Nielsen  and  the  one  of  Tucker  were  limited 
to  two-dimensional  flows,  but  the  former  was  studied  because 
of  its  potential  capabilities  of  boundary-layer  calculation 
in  separated  flows,  and  the  latter  because  of  its  reliability. 
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Figure  4.  Comparison  of  the  viscous  methods. 
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SO  as  to  provide  a data  base  for  comparison  with  the  other 
methods.  The  hypothetical  body  shape  consisted  of  a cylinder 
with  a conical  afterbody  and  a conical  plume  expansion,  with 
the  downstream  plume  shape  being  assumed  a cylinder. 

The  first  three  viscous  programs  (Refs.  29-31)  are 
integral  methods,  differing  in  the  formulation  of  the  ordinary 
differential  equations  and  the  assumptions  inherent  in  the 
modeling  of  the  turbulent  viscous  parameters.  The  Kuhn- 
Nielsen  Turbulent  Separated  Boundary  Layer  Program,  however, 
provides  an  extension  of  the  usual  solution  into  the  realm  of 
separated  flows  through  rearrangement  of  the  independent  and 
dependent  variables,  thus  eliminating  the  problems  encountered 
with  singularities  in  the  boundary  layer  equations  near 
separation.  The  program  was  not  used  in  this  "separated 
mode,”  however,  for  the  test  case.  The  methods  of  Mayne  and 
Whitfield  (32,33)  both  are  based  on  the  work  of  Patankar  and 
Spalding  (34),  but  they  differ  primarily  in  their  treatment 
of  the  turbulent  viscosity,  the  former  using  an  eddy  viscos- 
ity model  while  the  latter  using  a turbulent  kinetic  energy 
approach.  Recent  work  of  Chou  (35)  has  indicated  that  the 
effects  of  transverse  curvature  may  be  substantial  for  body 
fineness  ratios  of  the  magnitude  associated  with  the  present 
configurations  considered,  so  the  Mayne  program  calculations 
were  perforned  not  only  with  the  two-dimensional  equations, 
but  also  with  the  axisymmetric  equations  with  and  without 
transverse  curvature. 
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As  can  be  seen  in  Figure  4,  page  27,  all  of  the 
methods  agree  well  in  the  prediction  of  the  displacement 
thickness  distribution  when  the  two-dimensional  equations  are 
used  (the  results  of  the  Whitfield  calculation  are  not  pre- 
sented as  they  agreed  very  well  with  the  other  Patankar  and 
Spalding  calculation).  With  the  axisymmetric  formulation, 
the  results  of  Bartz  and  Mayne  agree  fairly  well.  In  the 
Mayne  calculation,  the  effect  of  transverse  curvature  lessens 
the  axisymmetric  displacement  thickness  rise  near  the  after- 
body/plume junction,  but  compared  to  the  difference  between 
the  axisymmetric  and  two-dimensional  predictions,  the  differ- 
ence was  not  thought  to  be  especially  significant. 

As  a result  of  this  test  case  it  was  decided  that 
the  Kuhn-Nielsen  and  Tucker  program,  though  both  were  accu- 
rate, could  not  be  considered  because  of  their  two-dimensional 
limitations.  Further  investigations  with  the  remaining  programs 
involving  substantial  adverse  pressure  gradients  near  the 
afterbody/plume  junction,  which  were  realistic  physically,  indi- 
cated that  the  integral  method,  the  Bartz  program,  was  more 
capable  of  overcoming  these  gradients  than  either  of  the  Patankar 
and  Spalding  finite-difference  methods.  Since  extremely  steep 
pressure  gradients  can  occur  in  the  initial  iterations  of  the 
viscous/inviscid  method,  it  was  decided  that  this  characteristic 
of  the  Bartz  program  demanded  its  use  in  the  combined  viscous/ 
inviscid  formulation.  No  direct  comparisons  of  the  displacement 
thickness  were  made,  however,  between  the  Bartz  and  the^ Mayne 
programs  under  adverse  pressure  gradients  for  the  purpose  of 
judging  the  accuracy  of  the  programs. 
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The  Bartz  solution  is  based  on  the  numerical  inte- 
gration of  the  momentum  and  energy  equations  in  the  forms 


r r I \2''1/2 

11  = ^ 

^ J L 


drl^"'^^^  J 2-M^  + (aye)  dM  . l dr 

3?)  2 ^ ^ 

J L M 1 + 


respectively.  The  relationship  for  the  skin  friction  coeffi 
cient,  Cj,  is  based  on  a correlation  developed  by  Coles,  and 
is  derived  in  Appendix  A,  Bartz  (31) : 


C 


f 


; T 1 -0.6 

0.0256  1 aw 


The  Stanton  number,  C^,  is  represented  by  the  relationship 
based  on  the  Prandtl  number  correction  to  the  Reynolds 
analogy  due  to  von  Karman,  which  is  also  derived  in  Appendix 
A of  Reference  (31),  of  the  form 
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Once  these  two  relationships  are  assumed,  the  inte- 
gration of  the  momentum  and  energy  equations  may  proceed  once 
the  Mach  number  distribution  at  the  edge  of  the  boundary 
layer  and  the  initial  shape  parameter,  6*/6,  are  specified. 
The  former  distribution  is  provided  by  the  inviscid  calcula- 
tion, and  the  latter  by  assuming  a velocity  and  stagnation 
temperature  distribution  through  the  boundary  layer.  Bartz 
uses  a l/7th-power  law  distribution  for  both  the  velocity  and 
stagnation  temperature  variation  (actually,  - T^) . 

Typical  run  time  for  the  Bartz  program  along  the 
entire  length  of  a body  under  a varying  pressure  distribution 
is  about  twenty-five  CPU  seconds  on  the  IBM  370  mod  165. 

Plume  Calculation 

In  order  to  analyze  the  flow  over  a boattail  with 
a real  exhaust  plume,  it  is  necessary  to  be  able  to  determine 
the  "shape"  of  the  plume,  recognizing  that  the  shape  must  be 
a simplification  of  the  quite  complex  flow  actually  encoun- 
tered (36) . When  viscous  effects  are  considered,  the  actual 
shape  of  the  plume  is  open  to  question,  and  one  may  choose  to 
define  the  shape  using  a "dividing-streamline"  measure  or  a 
certain  percentage  of  the  difference  between  the  plume  and 
external  velocities,  to  cite  two  examples.  If  viscous 
effects  are  ignored,  however,  the  merging  of  the  plume  and 
external  flows,  both  considered  to  be  inviscid,  determine  a 
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well-defined  boundary  between  the  two  flows.  This  boundary 
is  referred  to  as  the  "inviscid  plume  boundary"  or  "inviscid 
reference  line." 

The  importance  of  the  inviscid  plume  boundary  is 
that  (1)  it  serves  as  the  basic  solid  body  shape  of  the  plume, 
which  is  the  basis  for  the  aforementioned  solid-body  or  dis- 
placement effect  of  the  plume,  and  (2)  it  is  the  reference 
line  to  which  corrections  to  the  plume  shape  due  to  entrain- 
ment effects  are  applied. 

The  plume  shape  is  calculated  using  a Method  of 
Characteristics  (MOC)  digital  computer  program  produced  by 
the  Lockheed  Missiles  and  Space  Company  (37) . This  program 
is  one  of  a group  of  programs  used  by  Lockheed  to  analyze 
various  aspects  of  the  nozzle  exhaust  plume  flow  problems  as 
well  as  external  transonic  flows  and  forces  and  moments  on 
objects  immersed  in  such  flows.  The  Lockheed  MCX:  program  used 
in  these  studies  is  valid  for  ideal  or  real  gases  in  super- 
sonic compressible  flow,  and  two-dimensional  or  axisymmetric 
geometry  may  be  used.  It  has  a number  of  starting  options  for 
the  flows,  and  either  solid  walls  or  free  boundaries  may  be 
specified  downstream  of  the  starting  point,  with  capability 
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of  switching  from  a solid  boundary  to  a free  boundary  at  the 
user's  option.  This  characteristic  makes  the  program  ideal 
for  those  cases  where  it  is  desirable  to  begin  the  plume  cal- 
culation within  the  nozzle,  for  example  at  the  sonic  line  in 
the  throat  of  a converging-diverging  nozzle,  and  to  continue 
it  past  the  lip  at  the  exit  of  the  nozzle.  In  cases  where 
the  nozzle  is  sonic,  that  is,  the  throat  of  the  nozzle  is  at 
the  exit  plane,  the  calculation  may  be  started  at  this  point 
and  the  flow  is  calculated  immediately  into  the  free-boundary 
region. 

When  used  to  furnish  the  plume  shape  for  the 
viscous/inviscid  iteration,  the  program  is  supplied  the  nozzle 
pressure  ratio  of  the  exhaust  and  external  flows,  i.e.  the 
nozzle  total  pressure  divided  by  the  free-stream  static  pres- 
sure, the  exhaust  total  temperature,  pertinent  gas  constants, 
the  Mach  number  distribution  of  the  starting  line,  and  the 
proper  boundary  conditions.  The  calculation  runs  less  than 
ten  seconds  CPU  time  on  the  IBM  370  mod  165.  Some  examples 
of  actual  typical  plume  shapes  are  shown  in  Figure  5.  The 
solid  lines  represent  the  plume  as  calculated  by  the  Lockheed 
program,  but  the  plumes  were  extended  as  shown  by  the  dashed 
lines.  The  difference  in  distribution  on  the  boattail 
caused  by  this  simplification  of  plume  shape  was  minor. 
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Combined  viscous/Inviscid.  Iterative  Program 

The  general  method  for  the  iterative  solution  of 
the  type  previously  described  is  shown  in  the  flow  chart  in 
Figure  6,  As  can  be  noticed  immediately,  the  technique  in- 
volves two  main  interation  loops,  an  "inviscid/inviscid" 

(I/I)  loop  and  the  "viscous/inviscid"  (V/I)  loop.  The  former 
loop  is  an  iteration  between  the  external  and  exhaust  inviscid 
flows  with  the  purpose  of  modifying  the  inviscid  plume  shape 
according  to  the  pressure  distribution  on  its  surface  as 
calculated  by  the  external  inviscid  program.  When  the  plume 
shape,  which  is,  in  the  ideal  case,  calculated  from  a longi- 
tudinal pressure  distribution,  causes  the  external  inviscid 
program  to  return  essentially  the  same  pressure  distribution, 
the  iteration  is  considered  to  be  converged.  At  this  point 
the  V/I  loop,  which  has  been  described,  is  entered.  Upon  con- 
vergence of  the  V/I  loop  further  refinements  of  the  solution 
may  be  carried  out  by  returning  to  the  I/I  loop  with  the  new 
effective  body  shape  calculated  in  the  V/I  loop. 

The  Combined  Viscous  Inviscid  (CVI)  program  performs 
only  the  calculations  in  the  V/I  loop,  that  is  to  say,  the 
running  of  the  two  inviscid  programs  in  the  I/I  loop  and  the 
transferral  of  converged  data  to  the  V/I  loop  are  carried 
out  manually.  Manual  operation  of  this  portion  of  the  com- 
plete problem  is  not  an  undue  hardship,  for  in  the  many  cases 
where  the  body  to  be  analyzed  has  a solid  sting,  the  I/I  loop 
is  not  required,  and  in  real  plume  cases  it  has  been  found 
that  more  than  one  iteration  is  usually  not  necessary  for  a 
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Figure  6.  General  iteration  flow  plan. 
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realistic  converged  plume  shape.  In  fact,  the  plume  shape  as 
determined  from  the  nozzle  pressure  ratio  pertaining  to  the 
ambient  external  flow  is  usually  very  close  to  the  plume 
shape  developed  to  convergence  in  the  I/I  loop,  and  it  is 
used  in  the  calculations  herein.  The  following  discussion 
will  therefore  concentrate  on  input  requirements  for  a given 
problem  and  the  operation  of  the  CVI  program,  and  the  con- 
verged invisid  plume  shape  (or  solid  sting  shape)  will  be 
assumed  known. 

When  confronted  with  a new  body  configuration,  the 
South-Jameson  Transonic  Program  (RA.XBOD)  is  usually  run  a 
number  of  times  for  the  purpose  of  obtaining  an  inviscid 
solution  for  the  body  with  no  boundary  layer.  RAXBOD  has  a 
number  of  parameters  which  control  the  operation  of  the  pro- 
gram, and  it  is  thus  advisable  to  "tune"  the  program  itself 
before  embarking  on  a full  run  with  the  CVI  program.  Table 
IV  lists  the  most  important  input  requirements  for  RAXBOD 
along  with  typical  values.  Parameters  that  control  output 
functions,  such  as  card  punching  or  plotting,  and  parameters 
that  remain  the  same  for  most  runs  are  not  included  in  this 
table . 

A brief  elaboration  on  some  of  the  parameters  pre- 
sented in  Table  IV  is  in  order  because  of  their  importance  in 
the  operation  of  RAXBOD.  The  stretching  parameters  at  the 
body  nose,  DSDXIO  and  BZ , are  adjusted  so  as  to  provide  a 
nearly  square  grid  at  that  location.  Failure  to  do  this  may 
result  in  unrealistic  solutions  in  the  area  of  the  nose  or 
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TABLE  IV 

INPUT  TO  THE  SOUTH-JAMESON  TRANSONIC  PROGRAM  (RAXBOD) 


Parameter 
N ame  ( s ) 

Explanation 

Typical 
Value (s) 

KKN 

Number  of  {XP,yP)  body  input 
coordinates 

97 

XP  ,YP 

(XP,YP)  body  coordinates  XP=0.0 

70.6, 

yp=o .0 
2.3,2 

71.4 

,0.1/«  * * r 

.3 

I MAX , JMAX 

Initial  number  of  XP  and  YP  grids 

25,25 

MIT 

Maximum  number  of  iterations 

200 

MHALF 

Number  of  mesh  halvings 

2 

COVERG 

Convergence  limit  of  transonic 
disturbance  potential 

0.001 

QF3 

Relaxation  parameter 

0.5 

DSDXIO 

Stretching  parcimeter  for  the  tan- 
gential coordinate  system  at  the 
body  nose 

1.0 

BZ 

Stretching  parameter  for  the  normal 
coordinate  system  at  the  body  nose 

0.0135 

DSDXIM 

Stretching  parameter  for  the  tan- 
gential coordinate  system  in  the 
afterbody  region 

1.0 

XI 

XP  location  of  beginning  of  sheared 
cylindrical  coordinate  system 

35.6 

SMAX 

Maximum  body  length  coordinate 

71.0 

GAM 

Ratio  of  specific  heats 

1.4 

AMINF 

Free-stream  Mach  number 

0.8 

^MIT  halved  at  each  mesh  halving 


^Convergence  occurs  when  - COVERG,  where 

is  defined  in  Table  II/  page  24 

^max 
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even  lack  of  convergence.  The  parameter  QF3,  which  is  the 
relaxation  parameter  in  the  South-Jameson  finite-difference 
scheme,  assists  convergence  when  it  is  close  to  or  equal  to 
0.5  at  a small  penalty  in  CPU  time.  Lower  values  of  QF3  may 
radically  effect  the  capability  of  converging  for  a given 
body  shape  and  Mach  number.  Finally,  the  sheared  cylindrical 
coordinate  system,  controlled  by  XI , should  be  begun  ahead  of 
any  concavity  in  the  body  shape  to  eliminate  any  possibility 
of  body-normal  coordinate  axes  intersecting  at  large  distances 
from  the  body. 

Once  an  acceptable  "body-alone”  solution  has  been 
obtained  from  the  RAXBOD  program,  the  CVI  program  may  be  used 
for  the  viscous/inviscid  iterative  calculations.  Required 
inputs  for  this  program  are  those  of  the  RAXBOD  program  plus 
additional  parameters  to  be  used  by  the  viscous  portion  of  the 
iteration.  The  most  important  of  these  parameters  are  pre- 
sented in  Table  V.  It  will  be  noted  that  there  is  some  dupli- 
cation of  parameters  between  the  inviscid  and  viscous  input, 
but  this  was  retained  so  that  if  difficulties  were  encountered 
in  the  CVI  program,  the  pertinent  input  data  could  be  trans- 
ferred directly  into  either  the  inviscid  or  viscous  program 
in  order  to  resolve  the  particular  problem. 

The  CVI  program  begins  its  iterative  procedure  with 
an  inviscid  calculation  of  the  flow  about  the  body  with  zero 
boundary  layer  displacement  thickness.  This  computation  is 
referred  to  as  the  "zeroth”  iteration,  as  it  is  a repeat  of 
the  inviscid  calculation  previously  made  using  RAXBOD  alone 
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TABLE  V 

INPUT  FOR  THE  VISCOUS  PORTION  OF  THE  CVI  PROGRAM 


Parameter 
Name (s) 

Explanation 

Typical 
Value (s) 

IXTAB 

Number  of  (X,Y)  body  input  co- 
ordinates 

97 

X.Y 

(X,Y)  body  coordinates 

X=0.0,0.2, . . . , 
70.6,71.4(in.) 
Y=0.0,0 .1, , . . , 
2.3,2.3(in.) 

MZETA 

Velocity  profile  power  law  expo- 
nent 

7 

ITVfTAB 

Wall  temperature  option: 

Adiabatic  = -1,  Constant  = 0, 
Table  = 1 

-1 

TO 

Free- stream  stagnationtemperature 

580  (‘>R) 

PO 

Free-stream  stagnation  pressure 

2000  (psf) 

GAMD 

Ratio  of  specific  heats 

1.4 

PRANDL 

Stagnation  Prandtl  number 

0.7 

ZMUO 

-5 

Viscosity  at  stagnation  temperature  1.29  x 10 

( Ibm/ft-sec) 

THETAI 

-4 

Initial  value  of  momentum  thickness  1.0  x 10 

(ft) 

PHII 

Initial  value  of  energy  thickness 

1.8  X 10"'^ 

(ft) 

EPSZ 

Geometry:  Asixymmetric  = 1, 

Two-dimensional  = 0 

1 
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during  the  "tuning"  operation.  The  purpose  of  the  zeroth 
iteration  is  to  provide  the  viscous  program  with  an  initial 
longitudinal  pressure  distribution  over  the  length  of  the 
body.  It  is  during  the  zeroth  iteration,  and  perhaps  the 
next  one  or  two  iterations,  that  unusually  steep  pressure 
gradients  and  extreme  pressures  may  develop,  and  the  boundary 
layer  displacement  thickness,  which  tends  to  ameliorate  the 
flow  expansion  and  compression  in  the  area  of  the  afterbody , 
is  either  non-existent  or  very  small. 


The  RAXBOD  program,  irrespective  of  the  (XP,YP)  body 
coordinate  pairs  furnished  as  input,  provides  the  converged 
pressure  distribution  at  longitudinal  locations  which  are 
functions  of  the  parameters  DSDXIO,  XI,  and  SMAX  (see  Table 
IV,  page  38) . Care  must  be  taken  that  the  {X,Y)  body  coordi- 
nate pairs  input  to  the  viscous  portion  of  the  CVI  program 
(Table  V)  correspond  exactly  to  the  RAXBOD  pairs,  which  have 


been  previously  determined  during  the  tuning  of  the  program. 
Slight  errors  in  the  longitudinal  location  of  the  pressures 
may  result  in  erroneous  converged  solutions  of  the  CVI  program, 
especially  in  regions  of  rapidly  changing  body  geometry.  Once 
the  boundary  layer  has  been  calculated  by  the  viscous  portion 
of  the  CVI  program,  the  (XP,YP)  coordinates  of  the  new  effec- 
tive body  are  transferred  back  to  the  inviscid  program  using 
the  following  relationships  (see  Figure  7) . 


XP  = 


Vd  + 


6*  sin  0 


bod 


“ \od 


6* 


^bod 
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and 


“ ’'bod  + «*  ®bod  = ’'bod  ^ «*  r 


where  the  longitudinal  derivative  of  the  body  shape  at  (X,Y) , 
that  is,  is  calculated  numerically.  Similarly,  the 

displacement  thickness  at  XP  is  then  calculated  by 


5*  = YP  - YRP, 

P 

where  the  body  radius  at  XP , i.e.  YRP,  is  calculated  by  inter- 
polation using  the  body  coordinates  that  were  input  to  the 
viscous  program. 

With  these  coordinates  of  the  new  effective  body, 
the  inviscid  calculation  is  repeated,  and  a new  longitudinal 
pressure  distribution  is  calculated.  This  point  is  then  the 
end  of  the  first  complete  viscous/inviscid  iteration,  and  the 
pressure  distribution  may  be  compared  with  that  resulting 
from  the  zeroth  iteration.  If  there  is  good  agreement  between 
the  two  pressure  distributions,  a highly  unlikely  occurrence, 
the  solution  is  considered  to  be  converged.  If,  as  is  more 
likely  the  case  at  this  early  stage  of  the  viscous/inviscid 
iteration,  there  is  little  agreement,  the  second  iteration  is 
begun  using  the  pressure  calculated  at  the  end  of  the  first 
iteration.  This  process  is  continued  until  convergence  of 
the  pressure  distributions  is  satisfactory. 

It  has  been  found  helpful,  and  in  some  cases  neces- 
sary, to  reduce  the  amount  of  the  displacement  thickness 
added  at  each  iteration,  thus  providing  a form  of  "relaxation" 
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of  the  incremental  displacement  thickness.  This  technique  is 
shown  schematically  in  Figure  8,  where  6*  represents  the 
displacement  thickness  that  was  used  in  the  (n-l)-st  viscous/ 
inviscid  iteration,  and  6*"  represents  the  displacement  thick- 
ness calculated  during  the  n— th  iteration  based  on  the  pres- 
sure distribution  from  the  (n-l)“st  iteration.  Rather  than 
use  the  whole  magnitude  of  ® "relaxed"  displacement 

thickness  is  defined: 

5*n  = + R - <5*"“^), 

p rel  p P P 

where  R is  a relaxation  factor, 

0 < R < 1, 

which  is  usually  between  0.3  and  0.6,  depending  on  the  sever- 
ity of  the  body  geometry  and  flow  conditions. 

The  necessity  for  such  a relaxation  arises  because 
of  the  sensitivity  of  the  transonic  flow  to  small  changes  in 
body  shape  and  the  resulting  effect  on  the  viscous  calcula- 
tions. 

An  extreme  pressure,  say  at  the  afterbody  plume  (or 
sting)  junction,  may  produce  a very  large  displacement  thick- 
ness from  the  viscous  program  calculations  which,  at  the  next 
inviscid  computation,  may  then  result  in  too  low  a pressure. 
This  low  pressure  may  then  cause  an  unrealistically  low  dis- 
placement thickness,  causing  too  high  a compression  in  the 
next  inviscid  calculation,  and  then  once  again,  too  large  a 
displacement  thickness  at  that  point.  The  solution  will  then 


44 


AEDC-TR-77-106 


AEDC-TR-77-106 


tend  not  to  converge,  but  to  oscillate  between  two  distinct 
solutions.  Relaxation  of  the  displacement  thickness  incre- 
ment will  greatly  aid  in  reducing  the  probabilities  of  such 
tendencies . 

Figure  9 shows  the  maximum  difference  in  along 
the  body  between  a given  iteration  and  the  previous  iteration 
for  four  values  of  relaxation  factor,  R.  The  body  being  con- 
sidered is  a sting-mounted  cylindrical  model  with  a circular 
arc  afterbody  of  length  to  maximum  diameter  ratio  of  1.768, 

(38) . The  convergence  of  the  solution  with  R = 1.0,  which 
corresponds  to  no  relaxation,  is  doubtful  and  the  oscillation 
between  two  solutions  is  implied  by  the  regular  sawtooth 
pattern  above  = 0 starting  with  the  thirteenth  itera- 

tion. A relaxation  factor  of  0.7  improves  the  situation,  but 
there  is  a lightly  damped  oscillation  up  to  the  fifteenth 
iteration.  For  values  of  R of  0.5  and  0.3,  the  solution  con- 
verges rapidly  and  definitely  with  no  tendency  to  oscillate 
in  a regular  pattern.  Both  of  these  relaxation  factors  would 
be  acceptable  from  the  standpoint  of  a relatively  quick  and 
reliable  solution,  but  the  R = 0.5  case  seems  to  be  slightly 
quicker  in  converging.  It  should  be  mentioned  that  the  three 
lower  values  of  R converge  to  the  same  solution. 

Figure  10  shows  the  pressure  coefficient,  the  effec- 
tive body  shape,  displacement  thickness  as  computed  by  the 
viscous  program,  and  the  build-up  of  the  relaxed  displacement 
thickness  in  the  area  of  the  afterbody  and  sting  by  a value 
of  R = 0.3.  At  convergence,  the  relaxed  displacement 
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thickness  attains  the  value  of  the  actual  displacement  thick- 
ness along  the  body,  and  this  may  be  seen  at  high  iteration 
numbers  in  Figure  11,  In  this  same  figure,  the  relaxation 
factor  may  be  obtained  from  the  slope  of  the  first  iteration 
line.  Figure  12  shows  the  magnitude  of  the  oscillating  pres- 
sure coefficient  and  relaxed  displacement  thickness  for  the 
case  of  R = 1,0  relative  to  the  well-converged  R = 0,3  case. 

It  can  be  seen  that  the  R = 1.0  solutions  for  the  nineteenth 
and  twentieth  iterations  bracket  the  converged  R = 0.3  solu- 
tion, and  convergence  of  the  R = 1.0  solutions  toward  this 
solution  is  extremely  slow,  if  at  all. 

The  criterion  for  convergence  of  the  viscous/invis- 
cid  iteration  is  usually  assumed  to  be  - 0.001,  as  this 

is  on  the  order  of  one-tenth  of  one  percent  of  the  stagnation 
pressure  coefficient  and  is  also  less  than  one  percent  of  the 
pressure  coefficient  range  in  the  area  of  the  afterbody  and 
plume.  However,  each  case  is  examined  individually  to  ensure 
that  the  solution  exhibits  no  tendencies  toward  divergence. 

The  CVl  Program  is  neither  a small  program  with 

regard  to  computer  core  storage  nor  a fast  program  with  regard 

to  CPU  run  times.  By  far  the  largest  portion  of  the  core 

storage  and  run  time  requirements  are  attributable  to  the 

inviscid  portion  of  the  iteration,  RAXBOD.  It  has  been  found 

that  the  tight  convergence  requirement  quoted  by  South  and 

Jameson  on  A4  (see  Table  II,  page  24)  may  be  reduced  dur- 
rri3-x 

ing  the  course  of  the  viscous/inviscid  iteration  from  0.0001 
to  0.01  with  no  appreciable  change  in  the  converged  solution. 
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Figure  11.  Boundary  layer  relaxation  displacement  thickness 
development  relative  to  calculated  displacement  thickness. 
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Figure  12.  Comparison  of  the  oscillating  solution  of  an  R = 1.0 
case  and  the  converged  solution  of  an  R = 0.3  case. 
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Even  with  this  and  similar  methods  of  making  the  CVI  Program 
more  efficient,  the  version  of  the  program  currently  in  use 
requires  540  kilo-bytes  of  storage,  and  a typical  case  will 
use  approximately  30  to  45  minutes  of  CPU  time. 
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CHAPTER  III 

MIXING  ANALYSIS  AND  THE  PLUME  ENTRAINMENT  MODEL 
Introductory  Remarks 

As  was  described  in  the  beginning  of  Chapter  II , 
the  effect  of  a real  plume  on  the  external  flow  over  an  after- 
body may  be  considered  as  consisting  of  two  components: 

(1)  a blockage  or  displacement  effect,  which  causes  the  exter- 
nal flow  to  react  to  the  plume  by  deflecting  as  if  it  were  a 
solid  body,  and  (2)  an  entrainment  effect,  which  accelerates 
the  external  flow  through  the  turbulent  momentum  exchange 
mechanism  between  the  two  flows.  If  the  first  effect  is 
treated  by  considering  the  plume  to  be  a solid  body,  as  shown 
in  Figure  2,  page  14,  then  the  second  effect  may  be  approxi- 
mated by  a suitable  modification  of  this  solid  body  plume 
shape.  It  is  the  purpose  of  this  chapter  to  derive  the  neces- 
sary relationship  for  the  calculation  of  an  effective  plume 
boundary  layer  displacement  thickness  correction,  6*,  relative 
to  the  known  Inviscid  Reference  Line  (IRL) . 

The  method  of  analysis  used  to  calculate  the  dis- 
placement thickness  correction  is  based  on  the  mixing  analysis 
usually  attributed  to  Chapman  and  Korst.  The  development  of 
this  method,  in  which  emphasis  was  placed  on  the  base  flow 
problem,  is  presented  in  References  (39-45).  Most  of  the 
early  efforts  assumed  no  initial  boundary  layer  profile  which 
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led  to  fully  developed  solutions  with  velocity  profiles  of 
the  error  function  type.  Hill  (46),  following  the  general 
method  of  Chapman  and  Korst,  took  the  initial  boundary  layer 
effect  into  account  in  the  base  flow  problem  in  two  ways.  In 
the  first,  he  assumed  fully  developed  velocity  profiles  for 
simplicity  of  computation,  but  with  the  origin  of  the  coordi- 
nate system  displaced  so  that  the  "equivalent  mixing  flow"  had 
the  same  mass  and  momentum  fluxes  as  the  initial  boundary 
layer  flows.  In  the  second,  he  assumed  a linear  variation  of 
the  eddy  viscosity  but  with  no  assumption  of  fully  developed 
velocity  profiles.  Both  methods  yielded  satisfactory  results. 
Bauer  and  Matz  (47)  approached  the  initial  boundary  layer 
problem  by  applying  the  principle  of  conservation  of  momentum 
below  the  dividing  streamlines  of  the  two  flows  and  Prandtl's 
Mixing  Length  Theory  to  obtain  a streamwise  distribution  of 
the  shear  stress.  Bauer's  Integral  Method  (48)  was  then  used 
for  the  solution  of  the  resulting  relationships.  Bauer  and 
FOX  (49)  have  extended  this  method,  using  the  turbulent 
kinetic  energy  (TKE)  method,  to  supersonic  nozzle  afterbody 
flows . 

The  method  used  in  the  following  analysis  is  based 
on  the  work  of  Childs  (50)  and  Korst  and  Chow  (51) , but  is 
extended  to  provide  the  effects  of  initial  boundary  layers  on 
all  pertinent  quantities . 

Mixing  Analysis 

The  basic  flow  model  to  which  this  mixing  analysis 
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was  applied  is  shown  in  Figure  13.  Two  flows,  both  with 
initial  boundary  layers,  begin  mixing  at  the  end  of  an  in- 
finitely thin  flat  plate,  X = 0.  This  point  represents  the 
nozzle/afterbody  and  plume  junction.  Further  downstream,  at 
X = , after  sufficient  momentum  exchange,  the  velocity  pro- 

file attains  a smooth,  fully  developed  form  with  very  little 
indication  of  the  presence  of  any  initial  boundary  layers. 

The  fully  developed  velocity  profile  is  characteristic  also 
of  the  family  of  profiles  formed  by  the  mixing  of  two  uniform 
streams,  i.e.  with  no  initial  boundary  layers,  and  it  was 
using  this  simplified  form  of  the  velocity  distribution  that 
Korst  and  Chow  developed  their  analytical  results.  It  was 
recognized,  however,  that  the  initial  boundary  layers,  parti- 
cularly on  the  afterbody,  could  be  comparable  in  size  to  the 
dimensions  of  the  afterbody  and  thus  could  not  be  ignored. 

This  analysis  will  retain  the  initial  boundary  layer  terms  in 
the  velocity  profile  function  and  other  quantities  developed 
in  the  derivation. 

The  most  serious  restriction  within  this  analysis 
is  that  of  two-dimensionality  of  the  flow.  It  was  determined, 
by  reviewing  the  results  of  previously  referenced  work  on  free 
turbulent  mixing,  that  the  displacement  thickness  correction 
applied  to  the  plume  would  be  small  in  comparison  with  the 
radius  of  the  plume  for  the  geometric,  aerodynamic,  and 
thermodynamic  parameters  typical  of  nozzle/afterbody  flows. 
Thus,  the  axisyrametric  effects  would  normally  be  of  a lower 
order  and  could  be  ignored  in  the  present  analysis. 
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Furthermore,  it  was  felt  that  the  availability  of  theoretical 
and  experimental  results  for  the  case  of  two-dimensional  mix- 
ing would  provide  a valuable  means  of  checking  the  present 
computation. 

The  analysis  is  valid  for  compressible,  turbulent, 
non-isoenergetic  attached  flow  with  a turbulent  Prandtl  number 
assumed  to  be  unity.  The  mixing  is  furthermore  assumed  to  be 
isobaric.  This  is  consistent  with  the  approach  of  Korst  and 
Chow,  but  is  not  fully  realistic  considering  the  longitudinal 
pressure  distributions  occurring  in  nozzle/afterbody  plume 
flows.  However,  even  though  the  pressure  at  the  nozzle/after- 
body and  plume  junction  may  become  quite  large,  recovery  to 
free-stream  conditions  is  usually  rapid.  It  was  concluded 
that  this  tendency  toward  rapid  pressure  recovery  and  the 
dominance  of  the  turbulent  mixing  mechanism  would  make  the 
isobaric  assumption  acceptable  or,  at  any  rate,  introduce  no 
more  error  than  the  constant  pressure  assumption  used  in  the 
calculation  of  the  inviscid  plume  shape. 

The  physical  coordinate  system,  shown  in  Figure  13, 
consists  of  the  orthogonal  (X,Y)  axes,  and  the  inviscid 
reference  line  (IRL) , or  plume  boundary,  is  a straight  line 
coincident  with  the  X-axis.  The  physical  coordinate  system 
is  aligned  with  the  reference  flow,  the  lower  flow  in  Figure 
13,  which  corresponds  to  the  nozzle  jet  flow.  A constant 
transverse  velocity,  v^,  is  assumed  in  the  upper  flow  at  a 
large  but  finite  distance  from  the  X-axis , corresponding  to 
mass  entrainment  of  the  external  flow.  Mo  corresponding 
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transverse  velocity  is  assumed  for  the  lower  flow,  as  the 
physical  axes  are  aligned  with  this  flow  and,  furthermore,  a 
large  and  finite  distance  in  this  direction  could  be  made  to 
lie  on  the  actual  nozzle  centerline,  about  which  the  flow  is 
symmetric.  The  transverse  velocity  is  required  to  provide 
sufficient  mass  to  attain  the  fully-developed  velocity  profile 
without  suffering  a velocity  deficit  elsewhere,  analogous  to 

the  results  of  Moulden  (52) . 

The  equation  of  motion  for  the  isobaric  mixing  may 

be  written  as 


u 


(5) 


where  s is  the  eddy  viscosity  for  the  turbulent  flow.  Dimen- 
sionless variables  are  introduced. 


u ,1  - X 


and  the  eddy  viscosity  is  assumed  to  be  of  the  form 


E = f('l') 


where  fCij*)  is  an  undetermined  function  that  approaches  unity 
as  and  e is  defined,  according  to  Gortler  (53)  , as 

T'  ' ^ 


= 


4o 


X (u^  + Uj^) 


If,  furthermore,  it  is  assumed  that  v is  negligible 
and  that  the  non-linear  coefficient  u can  be  approximated  by 
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U + U, 

a b 


which  is  consistent  with  the  small  perturbation  analysis  of 
Pai  (54),  then  Equation  (5)  becomes  the  heat  conduction  equa- 
tion with  a diffusivity  of  one: 


_ 3 <) 


where  ^ 


2a  J 


if 


The  boundary  conditions  applicable  to  this  equation 
are,  for  the  initial  velocity  profiles, 


for  -c»  < j;  ^ - 


'Mo,e)  = 


“b 


if(0,c)  = 


for  0 i i;  < 1 


9(0rC)  = 1 


for  1 < r.  < 


and  for  the  upper  and  lower  boundary  conditions 
ct)(C,C)  -*■  for  c + 

9 ( 5 , ^ ) * 1 for  c -<-  «> 


The  functions  initial 

boundary  layer  profiles  of  the  jet  and  external  flows,  respec- 
tively, and  are  represented  in  our  case  by  the  power-law 
velocity  profiles 
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*23*^'  ■ '*2b*^’  ■ ♦bl  ■ ij*- 


The  solution  of  Equation  (6)  subject  to  the  bound- 
ary conditions  of  Equations  (7)  is 


4>(n,n  ) = -j 


+ erf  (n  - n ) - erf  n + ^ 


n + 

- a 


2a 


n-n. 


nzi L-^  dB  t ^ 

I /F 


2b  h 


dB 


(8) 


1 a a 

where  ti  - -*■  

2»'r 


t)  = cn  ->  0 ^ 
^ p X 


for  large  values  of  x 


and  the  error  function  has  the  usual  definition 

n 

' 2 

erf  n = — ^ I dB 

/ii~  J 


The  velocity  profile  is  thus  represented  by  a family  of  curves 
in  n with  one  parameter,  pp-  The  first  three  terms,  within 
the  brackets,  represent  the  contribution  to  the  velocity  pro- 
file of  the  fully-developed  flow.  These  terms  predominate  at 
small  n (far  downstream).  The  last  two  terms,  containing 
the  integrals,  represent  the  contributions  of  the  lower  and 
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upper  initial  boundary  layers,  respectively.  These  terms  are 
most  important  at  large  rip  {near  .the  nozzle  exit)  , and  become 
negligible  far  downstream.  The  relationship  between  the 
terms  is  shown  for  an  intermediate  value  of  n in  Figure  14. 

r 

A quick  approximate  method  for  obtaining  the  integrals  is 
presented  in  Appendix  A. 

This  family  of  velocity  profiles  was  interpreted  by 
Chapman,  Korst,  and  subsequent  investigators  to  hold  in  an 
orthogonal  "intrinsic"  coordinate  system,  not  usually  coinci- 
dent with  the  physical  coordinate  system,  as  a result  of  the 
approximations  inherent  in  the  development  of  the  equation  of 
motion.  This  system  of  coordinates,  designated  (x,y)  in 
Figure  13,  page  56,  is  related  to  the  physical  coordinates  by 
the  relationships 

X ^ X 

y = Y + y^(x)  , 

where  yjj,{0)  = 0,  and  yj^(x)  may  be  determined  with  the  aid  of 
integral  relations  for  the  conservation  of  mass  and  momentum 
in  the  following  forms; 


o a 

1 pu  dY  + j pu  dY  + 

■^b  o 

X 

r 

* J Pb  V*  = 

J pu  dy  + y^  (p^u^  - PpUj^) 

o 

^Rb 
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and 


□ a 


-6, 


"^b  / Pb^b^^^  = / ^ ^m  ‘Pa'^a  ' 


Pb^J 


Multiplying  the  first  equation  by  u^  and  subtracting  from  the 
second  gives 


n„  = n 


m 'Ra  1-fi 


ii-cl) 


/ 


<> 


A-C‘(}) 

a 


7~I 


Rb 


'^Ra 

/ 

HRb 


<t> 


TT 

A-cy 


\ 


dn 


, „ , '*b  “b  ,2 


where  Crocco's  energy  relationship  for  the  stagnation  tempera- 
ture profile  has  been  assumed,  i.e. 


= ^ ^ + 
Ja 


<t)  - <j). 


and  the  momentum  and  displacement  thicknesses  of  the  flow  are 
known.  Using  the  auxiliary  integrals,  the  equation  for  rij^^ 
may  be  rewritten: 
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% = ■ iq 


I2 


+ Ti. 


Pb  .2  . ®a  . 

(ft,  + -K + 


^ ^ ^ ^ 


(9) 


An  expression  for  the  location  of  the  dividing  streamline 
(DSL) , which  defines  the  boundary  across  which  no  net  mass 
transfer  occurs,  may  be  developed  by  considering  the  mass 
flow  below  the  DSL; 


m 


pu  dY  + (Yr^ 


- 5J  = 


pu  dy 


After  some  manipulation,  similar  to  that  used  to  derive  the 

n relationship,  the  following  expression  involving  may  be 
m J 

written; 


^"^Ra^ 


I2  (HRa) 


- n. 


Pv  o 0V.  0=. 

,2  b . a 

-T  ’^hZZ  ^ -T 


(10) 


Leaving  Equations  (9)  and  (10)  for  the  moment,  it 
is  possible,  in  order  to  obtain  an  expression  for  an  effec- 
tive displacement  thickness,  6*,  due  to  plume  entrainment,  to 
consider  mass  conservation  above  the  dividing  streamline: 
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pudY  = 


j I dx 

pu  dy  - j b D 


J pu  dy  - J 


^Rh'^y 


m 


If  Vj^  is  assumed  constant,  the  above  expression  can 


be  rewritten  as 


r 


o. 

To. 


\ 


- c: 


♦b 


V 


1 - c 


- Vm  ♦b  C 


J 


where 


If  a correction  to 
a displacement  thickness  due 

X 


I 


dx,  is  defined  as 


the  plume  shape,  that  is  to  say, 
solely  to  the  entrained  mass. 


*b 


f 


then  Equation  (11)  becomes 


6*  ^ 


ii(nj) 


•Ob 

“■oa 


(l-C^) 


1 


(12) 


Equation  (12)  is  an  explicit  expression  for  the 
displacement  thickness  due  to  plume  entrainment,  and  as  such 
is  applied  to  the  plume  boundary  at  each  viscous/inviscid 
iteration.  Equations  (9)  and  (10)  are  explicit  equations  for 
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n and  which  are  required  in  Equation  (12).  The 

m i ] 

magnitude  of  may  be  expected  to  be  that  of  the  plume 
adjustment  which  would  result  if  the  inviscid/inviscid  loop 
were  used.  A listing  of  the  computer  program  used  in  the 
solution  of  these  equations  is  presented  as  Appendix  B.  A 
description  of  this  program  and  a discussion  of  the  required 
input  are  presented  in  the  next  siib-chapter . 

Digital  Computer  Program 

A short  digital  computer  program,  DELPL,  has  been 
generated  for  the  purpose  of  solving  the  equations  developed 
in  the  previous  sub-chapter.  The  program  is  used  as  an  inte- 
gral portion  of  the  Combined  Viscous /I nvi sc id  Program,  supply- 
ing values  of  the  displacement  thickness  due  to  plume  entrain- 
ment at  longitudinal  locations  along  the  plume  as  specified 
by  the  calling  subroutines.  This  adjustment  of  the  plume 
shape  is  performed  at  each  iteration  according  to  the  re- 
lationship below: 

New  YP’’  = + 5* 

e 

where  the  value  of  YP  on  the  right-hand  side  is  the  value 
calculated  by  the  viscous  portion  of  the  CVI  Program  with 
the  plume  assumed  a solid  body.  Correspondingly,  the  dis- 
placement thickness  may  be  written  as 

New  6*^  = 6*"  + 6*. 
p p e 

6*  is  calculated  using  Equation  (12),  is  com- 

puted directly  from  Equation  (10),  and  Equation  (9)  is  used 
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for  ri^.  Notice  that  it  is  not  necessary  to  determine  n j , as 
only  the  value  of  the  integral  is  needed.  The  two 

auxiliary  integrals,  (no  ) and  Io(nTj  )r  are  calculated 
using  a 96-point  Gaussian  integration  formulation  (55)  . This 
quadrature  has  proven  to  be  of  very  good  accuracy,  being 
comparable  to  a trapzoidal  integration  of  one-thousand 
points,  as  determined  by  check  calculations. 

Computation  of  the  integrals  in  the  velocity  equa- 
tion is  performed  using  a closed- form  approximation  described 
in  Appendix  A.  This  approximation  provides  a decrease  in 
calculation  time  (CPU)  on  the  order  of  five,  with  little  or 
no  difference  from  these  integrals  when  they  were  evaluated 
using  the  96-point  Gaussian  integration. 

Parameters  input  to  DELPL  by  the  calling  subroutine 
are  shown  in  Table  VI.  A listing  of  the  program  is  given  as 
Appendix  B.  The  calculation  of  6*  in  DELPL  is  rapid  in  spite 
of  the  numerical  integrations,  each  calculation  of  6*  requir- 
ing  about  one  second  of  CPU  time. 

Stagnation  Temperature  Distribution  Near  the  Nozzle/Afterbody 

It  was  assumed,  in  the  development  of  the  equations, 
that  the  turbulent  Prandtl  number,  PJ^^,  was  everywhere  equal 
to  unity.  Accordingly,  the  variation  of  stagnation  tempera- 
ture throughout  the  mixing  flow  field  was  assumed  to  be  linear 
with  respect  to  the  local  velocity,  following  the  relation- 
ship developed  by  Crocco  (56) , which  is  valid  for  turbulent 
flow  with  unity  turbulent  Prandtl  number.  In  particular,  the 
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TABLE  VI 

INPUT  FOR  THE  PLUME  DISPLACEMENT  THICKNESS  SUB- 
PROGRAM (DELPL) , AS  PROVIDED  BY  THE  CALLING  SUBROUTINE 


Parameter 
Name (s) 

Explanation 

Typical 
Value (s) 

XP 

Axial  locations  along  plume  X = 

78 

76.209, 

.163,  ..  . (in.) 

FB 

Velocity  ratio  of  the  two 
streams 

0.6 

TTA 

Stagnation  temperature  of  the 
reference  stream 

575. 0®R 

TTB 

Stagnation  temperature  of  the 
secondary  stream 

575. 0°R 

CPA 

Gas  constant,  C , of  the  reference 
stream  ^ 

6006.0 

UA 

Velocity  of  the  reference  stream 

1175.0 

ft/sec 

DSTRA 

Displacement  thickness  of  reference 
stream 

0.0125  ft 

DSTRB 

Displacement  thickness  of  secondary 
stream 

0.0250  ft 

EN 

Velocity  profile  power  exponent 

7 

N 

Number  of  points  to  be  calculated 
(100  maximum) 

10 

rwR 

Output  control  parameter 

1 

DSTRE 

Values  to  the  plume  displacement 
thickness  returned  from  DELPL 

DSTRE  = 0, 
-0-008, 
-0.017,  - . . 

(in.) 
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relationship  chosen  by  Korst  and  Chow  (51)  was 


which  provided  the  linear  variation  with  appropriate  values 
of  the  two  free-stream  stagnation  temperatures  at  large 
lateral  distances  from  the  dividing  streamline.  This  linear 
relationship,  while  valid  in  the  flow  far  downstream,  becomes 
inappropriate  near  the  nozzle/afterbody,  where  the  stagnation 
temperature  characteristics  are  more  those  of  boundary  layer 
flows  than  free  mixing.  That  is  to  say,  one  would  expect  the 
stagnation  temperature  distribution  to  retain  its  boundary 
layer  characteristics  for  some  distance  downstream,  gradually 
changing  to  the  Crocco  linear  form  as  the  flow  becomes  more 
and  more  fully  developed.  The  difference  is  shown  qualita- 
tively in  Figure  15,  where  the  sketch  of  the  stagnation 
temperature  distributions  of  the  boundary  layer  flows  over  an 
adiabatic  wall  is  taken  from  Shapiro  (57)  . The  stagnation 
temperature  distributions  are  characteristic  in  that  they 
show  an  "overshoot"  near  the  boundary  layer  edge,  and  that 
they  attain  the  adiabatic  wall  temperatures  at  the  wall,  which 
are  not  necessarily  the  values  obtained  from  Equation  (13)  at 
(j)  = 0.  These  phenomena  were  analyzed  by  van  Driest  (58)  , who 
assumed  a variation  in  the  turbulent  Prandtl  number  across 
the  boundary  layer,  subdividing  the  boundary  layer  into  three 
regions.  A concise  summary  of  a number  of  analytical  studies 
devoted  to  this  problem  is  provided  by  Whitfield  (33) . 
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The  inadequacy  of  the  Crocco  relationship  became 
apparent  during  the  initial  calculations  of  the  mixing  pro- 
gram, DELPL.  Far  downstream,  the  results  of  the  numerical 
solutions  of  the  integrals  involved  in  the  calculations  and, 
indeed,  the  effective  displacement  thickness,  S*,  agreed  well 
with  the  computations  of  Korst  and  Chow  (51)  for  fully 
developed  flow.  Close  to  the  beginning  of  the  mixing  region, 
however,  there  was  a small  but  persistent  discrepancy  in  the 
value  of  these  quantities.  All  possible  sources  of  error, 
such  as  truncation  errors  in  the  integration  schemes,  were 
eliminated  as  sources,  and  the  sole  remaining  source  was  the 
stagnation  temperature  ratio.  A,  which  appears  in  the  inte- 
grals ^nd  ^2^^Ra^‘  order  to  demonstrate  this,  it 

is  first  necessary  to  calculate  the  values  of  (rin  ) » 
^2^’^Ra^'  ' ^m'  were  expected  at  the  nozzle/ 

afterbody  junction,  that  is,  at  n 

P 


Beginning  with  the  definition  of  , i.e. 


^l^^Ra^  T 


Rb 


’^Ra 


Ob 


,2 

'a 


- c: 


tl-C‘H 


dn 


(14) 


'Rb 


it  was  recognized  that  at  the  nozzle/afterbody  junction, 
n -*■  the  velocity  profiles  of  the  two  as  yet  unmixed  flows 
would  be  essentially  the  initial  boundary  layer  distributions. 
The  lower  flow  from  the  nozzle  would  consist  of  a power  law 
distribution  from  n = 0 to  n = hp#  and  the  remaining  flow, 
from  Hp  to  would  be  at  a constant  velocity,  4)  = 1. 
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Similarly,  the  upper  (external)  flow  would  be  a power  law 

6b 

velocity  distribution  from  ti  = 0 to  n — “ •j—  the  re- 

a ^ 

mainin^  flow  to  n = '^Rb  constant  velocity,  <)>  = 

Dividing  the  integral  in  Equation  (14)  into  the  four 
appropriate  integrals  gives 


np—  1 ^.^2^2  J 

■^oa  ^ b '^Rb 


'^°b  p2  .2 

To=  ~ '*’b 

Ol 


dn 


5 — dn  + 


V 

P 


^-^a  <i’2b 


*2a  . . 

5 dn  + 

A-ct 


2 

2a 


dn , 


(15) 


Using  the  relationship 


A ^2  .2 

T Oa  - C,  ♦ 
’■a  ■ 0 " 1 - C? 


which  results  from  the  expression  for  constant  enthalpy  and 
the  isobaric  assumption.  Equation  (15)  may  be  written,  after 
integrating  the  first  and  fourth  integrals  directly,  as 


lim 

Hp  ^ « 


I, (ho,)  = - s- 


l'‘'Ra'  " ■ "^b  T:  '’p 

a 


0 

r 


pu 

Pa^a 


dn 


b 

a ^ 


TT  ^ 


^ ""  " ■ "p’ 


(16) 
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where  the  definitions  of  <(>2^  and  (j)2jj  have  replaced  them  in 
the  remaining  integrals.  For  Hp  “f  the  variable  of  inte- 
gration, Ti,  becomes  the  physical  coordinate,  Y,  under  the 
relationship 


n = s 


and  if  the  integrals  and  their  limits  are  modified  accord- 
ingly, Equation  (16)  becomes 


“ iRa  - 


^ ^ “^b  3“ 


(17) 


where  the  definitions  of  6*  and  6^  were  used  to  eliminate  the 
integrals  in  Equation  (16) . 

A similar  procedure  for  the  second  auxiliary  inte- 
gral, 12^'^Ra^'  gives 


lim 


^2^'^Ra^  “ ’^Ra  “ ’^p 


- n 


P Pp 


(18) 


Thus,  for  !!  ->■“»  both  It  and  l2^^Ra^ 

evaluated  directly,  since  all  of  the  quantities  on  the  right- 
hand  side  of  Equations  (17)  and  (18)  are  known. 

Substitution  of  Equations  (17)  and  (IB)  into  the 
previously  derived  expression  (Equation  (9),  page  64)  for  the 
displacement  of  the  intrinsic  coordinate  system,  gives 
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lim 

ri 

n “ m 

P 


{19) 


This  identity  is  physically  satisfactory,  as  the 
intrinsic  coordinate  system  is  not  displaced  at  the  beginning 
of  the  mixing  region.  Similar  direct  substitution  of  Equa- 
tions (17>,  (18),  and  (19)  into  the  expressions  for  the 
integral  Ij^tn^),  Equation  (10),  page  64,  and  Equation 

(12),  page  65,  yields 


lim 

P 


il(nj) 


'b 


(20) 


^1  = 0 (21) 

n cm 

P 

Equation  (21)  indicates  that  at  the  beginning  of 
the  mixing  region  there  is,  as  yet,  no  correction  to  the 
external  boundary  layer  displacement  thiclcness , 6*,  due  to 
entrainment. 

The  cases  analyzed  herein  are  all  of  the  "cold  flow” 
type,  that  is,  the  exhaust  is  provided  by  pressurized  air  in 
a reservoir  at  essentially  ambient  temperature  rather  than  by 
a combustion  process.  Thus,  for  our  cases,  may  be 

assumed  to  a close  degree  of  approximation,  and  the  Crocco 
relationship  reduces  to  A = 1 for  all  41.  This  simplified 
stagnation  temperature  distribution  along  with  a qualitative 
representation  of  the  actual  boundary  layer  temperature  dis- 
tributions are  shown  in  Figure  16.  A linear  approximation  to 
the  actual  stagnation  temperature  distribution  used  for  com- 
putation purposes  is  shown  as  a dashed  line.  The  parameter  C 
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is  a constant  evaluated  by  iteration  at  the  initial  station , 
that  is,  at  the  nozzle/afterbody  junction,  according  to  the 
modified  relationship 


lim 

n -►  « ^1^'^Ra^ 
P 


■Ob 


- (Jlr 


(l-cp^ 


(A-f^C)  - 


T-7 


dn 


= n 


Ra 


- n. 


where,  for  our  cases,  Toj^/To^  = A = 1 and  = 1 at  np  ^ 

Thus,  the  integral  is  evaluated  exactly  numerically 

at  the  initial  X-location,  which  is  usually  set  at  a high 

4 

value  of  Hp,  approximately  10  . In  order  to  evaluate  f ^ , it 
is  assumed  that  the  mechanism  that  produces  the  change  in  the 
stagnation  temperature  distribution  is  the  difference  in 
actual  velocity  profile  from  that  of  the  error  function  dis- 
tribution in  the  fully  developed  case.  The  necessary  function, 

f , may  be  written 
c 


0 - 4 - 

act  erf 


4 


act 


^ ^act  ” ^erf 
o o 


4 


erf 


r 

where  4=1  4 the  subscript  "act”  denotes  the 


76 


AEDCTR-77-106 


integral  of  the  actual  velocity  distribution  at  a given  sta- 
tion, the  subscript  "erf"  denotes  the  integral  of  the  fully 
developed  velocity  distribution  at  that  station,  and  the 
subscript  "o"  denotes  those  integrals  evaluated  at  the 
initial  station.  Thus,  f varies  from  unity  at  the  initial 
station  to  zero  far  downstream,  where  = '^erf'  inte- 

gration limits  aiTid  are  found  by  determining  the  values 
of  Ti  where  the  velocities  are  equal  to  0.9999  of  their  free- 
streatn  values. 

This  modification  to  the  A function  is  satisfactory 
in  obtaining  the  desired  results,  that  is,  the  effective  plume 
displacement  thickness,  6*,  at  the  initial  station  becomes 
essentially  zero.  A two-constant  approach  was  also  tried, 
in  which  not  only  but  also  l2^^Ra^  were  matched  with 

their  limiting  values,  but  the  slight  change  in  results  did 
not  justify  the  added  complication. 

The  Mixing  Parameter,  a 

One  of  the  most  critical  parameters  in  use  in  this 
analysis  is  the  mixing  parameter,  c,  also  known  as  the  spread- 
ing, spread  rate,  or  similarity  parameter.  Introduced  in  the 
expression  for  the  turbulent  viscosity,  (see  page  5G) , the 
mixing  parameter  is  essentially  a measure  of  the  rate  of  ^ 
spreading  of  the  exhaust  jet,  with  high  values  of  c indicating 
a low  spreading  rate.  Gortler' s (53)  original  analysis 
assumed  that  the  velocity  profiles  downstream  were  simply 
functions  of  (y/x) , hence  the  designation  "similarity 
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parameter."  Analyzing  the  case  of  two-dimensional  mixing  of 
a single  stream  and  a quiescent  atmosphere,  he  determined  the 
incompressible  value  of  ff  to  be  13.5,  based  on  the  data  of 
Reichardt  (59) . Further  investigations  of  the  incompressible 
single  stream  case  by  Leipmann  and  Laufer  (60)  based  on  their 
data  indicated  that  a value  of  o = 11  best  suited  the  data, 
while  in  comparison  with  their  data  using  Tollmien's  theory 
(61),  they  chose  a = 12.  Using  Tollmien's  data,  Cordes  (62) 
determined  o to  be  11.95.  Generally,  a value  of  o = 12  is 
accepted  as  the  "round-number"  value  for  the  incompressible 
one-stream  case.  Other  investigations  have  shown  that  a in- 
creases with  increasing  Mach  number,  with  the  relationship 
developed  by  Korst  and  Chow  (51)  being  accepted  as  in  reason- 
able agreement  with  available  data,  i.e. 

o = 12  + 2.758  M.  (22a) 

More  recent  efforts  in  determining  the  values  of  a 
for  the  base  flow  problem  using  a variety  of  turbulence  models 
have  been  summarized  by  Greenwood  (63)  in  his  doctoral  thesis. 
Since  the  emphasis  in  his  summary  is  the  effect  of  compres- 
sibility at  Mach  numbers  up  to  eight  on  the  single-stream  con- 
figuration, the  results  have  limited  applicability  to  the 
present  study. 

For  the  two- stream  mixing  case,  that  is,  where  both 
streams  have  non-zero  velocities,  Weinstein's  experiments  (64) 
showed  that  a varied  greatly  with  the  initial  velocity  ratio 
of  the  two  streams.  The  most  commonly  accepted  measure  of  this 
effect  is  approximated  by  a relationship  developed  by  Korst  and 
Chow  (51),  i.e. 
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a 


II 


1 + 


(22b) 


In  choosing  acceptable  models  of  a for  use  in  the 
mixing  portion  of  the  CVI  Program,  this  author  was  influenced 
by  the  pragmatic  viewpoint  of  Gauthier  in  his  justification 
for  the  choice  of  turbulence  models  for  his  numerical  bound- 
ary layer  program,  PETULA: 

Mais  il  ne  faut  pas  perdre  de  vue  que  nous^ 
poursuivons  un  travail  d’ingenieur,  c'est-a-dire 
que  notre  objectif  n'est  pas  de  contribuer  a 
une  amelioration  de  nos  connaissances  theoriques 
sur  la  turbulence,  mais  de  fournir  une  methode 
pratique  de  calcul  qui  soit  un  compromis  raisonnable 
entre  precision  et  cout.  (65) 

(However,  it  is  necessary  not  to  lose  sight 
that  we  are  pursuing  the  work  of  an  engineer, 
that  is  to  say,  our  objective  is  not  to  con- 
tribute to  the  improvement  of  our  theoretical 
knowledge  of  turbulence,  but  to  furnish  a 
practical  method  of  calculation  that  is  a 
reasonable  compromise  between  precision  and  cost  - 
author’s  translation.) 

It  was  decided  to  use  commonly  accepted  values  of  o 
in  the  viscous/inviscid  iterations  in  order  to  determine 
their  validity  in  the  cases  analyzed.  Accordingly,  each  con- 
figuration was  run  at  least  three  times  in  the  CVI  Program 
using  the  following  values  of  a; 

(a)  One-stream  o from  Equation  (22a) , 

(b)  Two-stream  a from  Equation  (22b), 

(c)  At  least  one  value  of  o empirically  chosen 
from  the  results  of  (a)  and  (b) . 
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Both  of  the  values  of  a in  (a)  and  (b)  are  further 
modified  for  the  effect  of  the  initial  boundary  layer  by 
means  of  the  relationship 


o_ 

0. 


x/(6a  + fifa) 

37 


1/7 


which  was  the  result  of  a curve- fit  of  data  from  Chapman  and 
Korst  (66),  Figure  17.  The  choice  of  a in  the  third  case, 

(c) , was  based  on  the  results  of  (a)  and  (b)  and,  as  such, 
indicated  the  magnitude  of  the  deficiencies  in  the  values  of 


a predicted  under  the  multitude  of  assumptions  and  simplfica 
tions  previously  described.  The  choices  for  (c)  will  be 
described  for  each  configuration  as  appropriate. 
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CHAPTER  IV 

ANALYTICAL  RESULTS  IN  COMPARISON  TO  EXPERIMENT 
Introductory  Remarks 

There  have  been,  over  the  course  of  the  years, 
numerous  wind  tunnel  tests  of  a variety  of  isolated  axisym- 
metric  body  shapes,  both  with  solid  stings  and  real  plumes, 
in  transonic  flows.  Many  of  the  configurations  would  be 
directly  applicable  for  comparison  with  the  flow  predictions 
of  the  CVI  Program,  Reference  to  the  recent  AGAEiD  efforts 
toward  improved  nozzle  testing  techniques  (67)  and  the  exten- 
sive bibliography  presented  by  the  U,  S.  Army  Missile  Command 
(68),  to  cite  just  two,  will  indicate  the  level  of  such  ef- 
forts in  recent  years. 

The  choice  of  configurations  to  be  analyzed  in  this 
chapter  was  strongly  influenced  by  the  works  presented  at  the 
Propulsion  Interactions  Workshop  at  Langley  Research  Center 
in  May  1976  (16)  and  recent  work  at  the  Arnold  Engineering 
Development  Center.  It  was  felt  that  such  configurations 
truly  represented  current  areas  of  interest  in  the  field, 
and  furthermore,  that  future  efforts  would  probably  be  direc- 
ted toward  similar  configurations  for  some  time  to  come.  This 
necessary  narrowing  of  the  field  left,  unfortunately,  a great 
number  of  works  untouched,  but  the  analyses  presented  may  be 
considered  representative  and  pertinent. 
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Solid  Sting  Cases 

Equivalent  body.  One  of  the  earliest  uses  of  the 
CVI  Program  was  in  the  analysis  of  a body  of  revolution  whose 
cross-sectional  area  distribution  was  designed  to  be  repre- 
sentative of  that  of  a typical  twin  jet  fighter  aircraft,  that 
is  to  say,  an  "equivalent  body."  The  wind  tunnel  tests  were 
carried  out  in  the  16-foot  Transonic  Wind  Tunnel  at  the  Arnold 
Engineering  Development  Center  (69)  , and  the  model  was 
equipped  with  interchangeable  contoured,  cylindrical,  and  15- 
degree  boat tails.  Only  the  contoured  boattail  was  analyzed, 
as  prime  interest  at  that  time  was  centered  about  the  capa- 
bility of  the  CVI  Program  to  predict  pressures  along  the 
entire  length  of  an  axisymmetric  body,  and  it  was  felt  that 
the  analysis  of  one  boattail  configuration  would  be  sufficient 
to  demonstrate  this  capability. 

As  Figure  18  shows,  the  equivalent  body  has  a some- 
what exotic  geometry,  particularly  in  the  areas  around  Sta, 

20,  Sta.  90,  and  along  the  boattail.  These  areas  first  pre- 
sented some  problems  for  the  inviscid  portion  of  the  CVI 
Program  as  the  body-normal  coordinates  tended  to  intersect  at 
great  distances  from  the  body,  but  these  problems  were 
remedied  by  moving  the  beginning  point  of  the  sheared  coordi- 
nate system  ahead  of  this  area.  Also  to  be  noted  is  the  step 
of  1.268  inches  at  the  end  of  the  contoured  boattail  (Sta. 
130.053).  Appealing  to  the  "wake-body"  concept,  in  which  an 
effective  body  shape  is  sought  in  order  to  represent  the 
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Figure  18.  Equivalent  body  shape  and  circumferential  pressure  tap 
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separated  region  after  the  step,  it  was  decided  to  modify 
this  area  by  means  of  a conical  surface  from  the  end  of  the 
boattail  to  the  sting.  The  angle  of  the  conical  surface  to 
the  horizontal  axis  was  chosen  to  be  7.5  degrees,  based  on 
extrapolation  of  the  Mach  number  dependent  analytic  results 
of  Mueller  (70)  and  the  transonic  results  of  Wu,  et  al.  (71) 
and  Chen  (72) . It  was  found  by  the  investigators  in  the  last 
two  references  that  this  angle  is  a weak  function  of  the  Mach 
and  Reynolds  numbers  in  transonic  flow.  This  angle  was  used 
for  all  Mach  numbers  analyzed  and  was,  in  fact,  little  more 
than  a realistic  estimate,  but  the  effect  on  the  flow  over 
the  equivalent  body  as  a whole  was  anticipated  to  be  negli- 
gible . 

The  equivalent  body  was  analyzed  at  three  Mach  num- 
bers: 0.598,  0.892,  and  1.196.  The  nominal  unit  Reynolds 

number  for  all  three  cases  was  2.5  x 10^  ft  The  results 
are  presented  in  Figures  19a-c.  In  the  two  subsonic  cases, 
Figures  19a  and  b,  the  agreement  between  the  predicted  and 
experimental  values  of  pressure  coefficient  are  extremely 
good.  All  of  the  variations  in  pressure  due  to  the  concave 
body  geometry  are  accounted  for.  The  CVI  Program  predicts 
too  high  of  a pressure  along  the  boattail,  in  the  area  of 
Sta.  123,  but  considering  the  rough  estimate  of  cone  angle 
used  in  the  wake-body  approximation,  agreement  is  still  very 
good.  In  the  supersonic  case,  Figure  19c,  it  is  apparent 
that  the  nose  region  is  not  well  predicted,  and  this  error 
results  in  deviations  from  the  measured  pressures  up  to  Sta. 
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a.  Mach  Number  = 0.598 

Figure  19.  Longitudinal  pressure  coefficient  distribution  for  the 
equivalent  body. 
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60.  Th£s  problem  can  sometimes  be  eliminated  by  changing 
the  grid  spacing  in  that  area.  Farther  downstream  the  agree- 
ment is  good,  and  the  CVI  Program  predicts  the  compression 
pressure  fairly  well,  lending  credence  to  the  wake-body  con- 
cept for  small  steps  in  supersonic  flow. 

NASA  Langley  circular  arc  boattail . Efforts  at 
NASA  Langley  Research  Center  in  recent  years  have  centered 
around  the  effect  of  jet  exhausts  on  boattail  pressure  drag. 
Pressure  and  drag  measurements  on  three  circular  arc  boat- 
tails  of  different  lengths  but  with  identical  closure  ratios 
of  0.50  were  carried  out  by  Reubush  and  Runckel  (73).  Pub- 
lished at  the  same  time  were  the  results  of  Reubush  (74)  on 
different  boattail  length  models  with  closure  ratios  of  0.6 
and  0.7.  Basic  conclusions  regarding  the  effects  of  boat- 
tail  length  and  closure  were  stated  in  these  reports, 

but  any  definite  correlations  were  obviated  by  the  very  strong 
effects  of  Mach  niamber  and  flow  separation  on  the  pressure 
distribution  over  the  boattail. 

The  emphasis  at  NASA  Langley  were  then  directed 
toward  determining  the  applicability  of  solid  plume  simula- 
tors for  afterbody  wind  tunnel  tests.  These  simulators,  in 
actuality  cylindrical  stings,  would  have  the  advantage  of 
providing  the  support  for  the  model,  thus  relieving  any 
necessity  of  a large  and  disruptive  strut,  and  would  also 
eliminate  the  need  of  extensive  internal  plumbing  required  by 
the  jet  exhaust  models  at  a great  reduction  in  cost.  The 
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results  of  the  experimental  investigations  were  published  by 
Reubush  (75) , and  a much  fuller  treatment  was  later  presented 
also  by  Reubush  (38) . It  is  from  this  last  document  that 
data  were  obtained  which  were  used  for  comparison  with  the 
CVI  Program  predictions  at  various  Mach  numbers. 

The  wind  tunnel  tests  in  Reference  (38)  were  per- 
formed in  the  Langley  16-foot  transonic  tunnel,  which  has  a 
Mach  number  range,  in  continuous  operation,  from  0.20  to 
1.30.  It  is  important  to  note  that  this  wind  tunnel  is 
atmospheric,  thus  each  Mach  number  differs  slightly  in 
Reynolds  number  (Reynolds  number  effect  will  be  compared 
later  with  data  from  Reference  (17)).  These  small  changes  in 
Reynolds  number  were  accounted  for  in  the  CVI  Program  by 
adjusting  the  input  parameters  in  the  viscous  portion  of  the 
program  accordingly. 

The  model  was  strut  supported,  and  the  plume 
simulator  was  eleven  inches  in  length,  giving  a length-to- 
model  maximum  diameter  ratio  of  1.83.  The  model  support 
installation  is  shown  in  Figure  20.  Eight  circular  arc  boat- 
tails,  varying  in  afterbody  length  and  closure  ratio,  were 
tested.  Of  the  eight  boattails  available,  the  one  with  a 
fineness  ratio  (£,/d  ) of  1.768  and  a closure  ratio  (d  /d  ) of 

Til  0X  m 

0.50  was  chosen  for  analysis.  A sketch  of  this  boattail  is 
shown  in  Figure  21.  This  particular  configuration  was  chosen 
because  its  fineness  ratio  permitted  the  flow  to  be  attached 
at  most  Mach  numbers,  and  yet  its  low  closure  ratio  allowed 
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Figure  20.  Sketch  of  the  NASA  Langley  1/6.  = 1.768  circular  arc 

boattail  model  showing  internal  sting  Srrangement  for  support 
of  simulators. 
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58.608  < X < 69.608 

Note;  All  dimensions  in  inches  in  the  above  equations. 

Figure  21.  Shape  of  the  NASA  Langley  Jl/djj,  = 1.768  circular  arc 
boattail  and  cylindrical  plume  simulator. 
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ample  opportunity  for  axisymmetric  effects  to  be  present. 

Results  of  the  CVI  Program  predictions  for  Mach 
numbers  ranging  from  0.4  to  1.3  are  presented  in  Figures 
22a-g.  Because  the  Langley  16-foot  wind  tunnel  is  atmos- 
pheric, the  unit  Reynolds  number  varied  with  Mach  number  from 
a value  of  2,40  x 10^  ft"^  at  = 0.4  to  4.22  x 10^  ft*^  at 
= 1,3.  Flow  predictions  are  presented  only  for  the 
boattail/sting  areas  of  the  model,  as  the  data  were  only 
available  for  this  portion  of  the  model.  The  CVI  Program, 
however,  utilized  the  full  model  and  sting  in  both  its 
viscous  and  inviscid  portions  during  the  course  of  the  itera- 
tions. 

In  general,  the  predictions  agree  well  with  the 
data.  Some  general  statements  may  he  made  concerning  the 
characteristics  of  the  converged  iterative  solutions. 

Firstly,  in  the  subsonic  cases,  the  expansion  of  the  flow 
over  the  shoulder  of  the  boattail,  from  Sta.  45  to  approxi- 
mately Sta.  55,  is  underestimated  slightly.  As  the  Mach 
number  is  increased,  however,  this  tendency  becomes  less,  and 
the  shoulder  predictions  are  very  good  at  the  higher  subsonic 
Mach  numbers.  Secondly,  at  the  afterbody /sting  junction,  the 
CVI  Program  generally  overpredicts  the  peak  compression  pres- 
sure, although  data  were  not  available  at  the  exact  point  of 
the  junction,  so  it  is  difficult  to  estimate  exactly  the 
magnitude  of  the  error.  The  prediction  of  a higher  peak 
pressure  at  this  location  may  indicate  a small  separated  area 
there,  as  its  magnitude  increases  as  Mach  number  is  increased. 
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Figure  22.  Longitudinal  pressure  coefficient  distribution  for  the 
NASA  Langley  £/d  = 1.768  circular  arc  boattail  in  the  Langley 

16-ft  transonic  ^ind  tunnel. 
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It  should  be  noted  that  for  Mach  numbers  greater  than  0.90 
the  shoulder  area  is  supersonic,  and  the  rapid  compression 
downstream  of  this  area,  but  upstream  of  the  boattail/sting 
junction,  suggests  the  formation  of  a weak  shock  wave. 

In  the  supersonic  case,  Figure  22g,  the  steep 
gradient  of  the  compression  indicates  the  presence  of  a strong 
shock  wave.  The  CVI  Program  predicts  the  expansion  over  the 
shoulder  quite  well,  differing  in  shape  as  it  does  from  the 
expansions  in  the  subsonic  cases.  After  the  shock  wave, 
however,  the  prediction  is  less  satisfactory.  The  tendency 
of  the  data  to  assume  a slightly  curved  plateau  aft  of  the 
boattail/sting  junction  indicates  an  area  of  flow  separation 
much  larger  than  that  found  in  the  subsonic  cases.  This,  of 
course,  moves  the  flow  out  of  the  realm  of  applicability  of  the 
CVI  Program.  Addition  of  a transonic  separation  criterion, 
such  as  found  in  Reference  (14),  could  be  valuable  in  this 
and  similar  cases. 

It  was  mentioned  before  that  the  length  of  the 
cylindrical  plume  simulator  was  eleven  inches.  The  effect  of 
this  short  sting  is  reflected  in  Figures  22  by  the  continued 
decrease  of  the  last  four  data  points  to  pressures  even 
lower  than  free-stream,  a result  of  the  flow  separation  and 
expansion  at  the  end  of  the  sting.  The  wake-body  model  was 
applied  in  a manner  similar  to  its  application  with  the 
equivalent  body,  that  is  to  say,  a 7.S-deg  cone  was  assumed  at 
the  end  of  the  cylindrical  plume  cylinder,  and  the  improvement 
in  the  predictions  is  shown  as  dashed  lines  for  Mach  numbers 
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Of  0.4,  0.8,  and  0.94. 

As  a point  of  interest.  Figure  23  shows  the  results 
of  the  CVI  Program  when  the  boundary  layer  was  calculated  as 
two-dimensional  for  a Mach  number  of  0.8.  Although  this  pre- 
diction is  about  as  good  in  the  shoulder  region,  when  the 
displacement  thickness  is  small,  as  the  case  with  the  axi- 
symmetric  boundary  layer  (Figure  22c,  page  96) , it  suffers  at 
and  downstream  of  the  boattail/sting  junction.  Clearly,  the 
axisymmetric  boundary  is  necessary  when  adverse  pressure 
gradients  and  changing  body  geometry  demand  a rapid  increase 
in  boundary  layer  displacement  thickness. 

To  measure  the  effect  of  Reynolds  number  changes  on 
the  pressure  distribution  over  this  boattail,  data  from  more 
recent  work  by  Reubush  and  Putnam  (17)  were  utilized.  The 
data  were  taken  at  three  azimuths  about  the  body,  but  differ- 
ences were  small.  The  wind  tunnel  used  for  these  tests  was 
the  NASA  Langley  1/3-meter  cryogenic  tunnel,  the  test  section 
of  which  is  essentially  a smaller  version  of  the  NASA  Langley 
16- foot  wind  tunnel.  By  varying  tunnel  stagnation  pressure 
and  stagnation  temperature,  Reynolds  numbers  (based  on  the 

length  of  the  model  from  the  nose  to  the  beginning  of  the 

G 6 

boattail)  from  5 x 10  to  over  50  x 10  are  obtainable.  Six 

models  were  tested,  and  the  model  chosen  for  analysis  was  the 

same  a/d  = 1.768  circular  arc  boattail  chosen  previously  in 
m 

this  section.  It  was,  of  course,  much  smaller,  being  one- 
sixth  the  size  of  the  model  used  in  the  NASA  Langley  16 -foot 
tunnel.  The  model  had  a cylindrical  plume  simulator,  but 
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Figure  23.  Longitudinal  pressure  coefficient  distribution  for  the 
NASA  Langley  £/d  = 1.768  circular  arc  boattail  using  a two- 

dimensional  bounSary  layer  in  the  CVT  program. 
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it  was  not  instrumented,  and  hence  no  sting  data  are  avail- 
able. The  results  of  the  CVI  Program  in  comparison  with  the 
data  are  shown  in  Figures  24a-f  for  three  Reynolds  numbers 
at  both  of  two  Mach  numbers,  0,6  and  0.9.  Although  there  is 
some  scatter  of  the  data,  particularly  in  the  area  of  the 
expansion  over  the  boattail  shoulder,  the  CVI  Program  pre- 
dictions agree  well  with  the  data.  For  both  Mach  numbers, 
the  effect  of  increasing  Reynolds  number  is  to  decrease  the 
pressure  at  the  shoulder  expansion  and  to  increase  the  com- 
pression at  the  boattail/sting  junction.  This  is  a result  of 
the  decrease  of  boundary  layer  displacement  thickness  as 
Reynolds  number  is  increased.  The  effective  body  shape  for  a 
high  Reynolds  number  will  be  closer  to  the  actual  body  shape, 
and  thus  the  pressure  distribution  will  tend  to  resemble  the 
body-alone  inviscid  solution  more  so  than  a low  Reynolds  number 
case.  Both  the  data  and  the  CVI  Program  predictions  exhibit 
this  characteristic,  but  the  effect,  even  with  a change  in 
Reynolds  number  of  almost  a magnitude,  is  very  small.  The  one 
exception  is  the  expansion  pressure  in  Figure  24e,  which 
appears  to  be  slightly  too  low.  This  may  be  due  to  a slight 
difference  in  test  Mach  number  (about  1%  higher)  for  this  case 
versus  the  higher  Reynolds  case.  Figure  24 f. 

CBA  Boattail.  The  CVI  Program  was  also  applied  to 
a boattail  configuration  with  a cylindrical  sting  similarly 
analyzed  by  Chow,  Bober,  and  Anderson  (9) , and  this  configura- 
tion is  denoted  by  the  acronyn  "CBA  Boattail.”  The 
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Figure  24.  (Continued) 
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experimental  results  used  for  comparison  were  those  of 
Shrewsbury  (76) . The  results  of  the  CVI  Program  were 
anomalous  at  high  Mach  numbers,  and  as  such  are  presented 
as  Appendix  C. 

Real  Plume  Cases 

NASA  Langley  circular  arc  boattail.  The  model  used 
for  the  first  real  plume  analysis,  the  NASA  Langley  ^■/dj^  = 
1.768  circular  arc  boattail,  was  identical  in  external  shape 
to  the  solid  sting  model  presented  earlier  in  this  chapter, 
page  89,  less,  of  course,  the  plume  simulator.  It  was  strut 
mounted,  and  the  internal  plumbing  arrangements,  used  to 
furnish  a "cold"  plume,  are  shown  in  Figure  25. 

Reubush  and  Runckel  (73)  tested  a number  of  con- 
figurations at  Mach  numbers  from  0.4  to  1.3  and  nozzle  pres- 
sure ratios  (NPR's)  from  one  (jet-off)  to  about  six.  The 
analysis  and  data  presented  herein  are  for  the  l/d^  = 1.768 
circular  arc  boattail  at  a Mach  number  of  0.8  for  NPR's  of 

approximately  two  and  four. 

The  high  pressure  air  was  provided  at  essentially 
ambient  conditions,  and  the  wind  tunnel  used  for  the  tests 
was  the  Langley  16-foot  transonic  tunnel,  which  was  atmos- 
pheric. It  was  thus  judged  reasonable  to  assume  the  ratio  of 
stagnation  temperatures  of  the  external  and  internal  flows  to 
be  unity  for  the  purposes  of  the  plume  and  mixing  calcula- 
tions . 

The  plume  was  calculated  using  the  Lockheed  Method 
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Figure  25.  Sketch  of  air-powered  cone-cylinder  model  with  a 
typical  circular-arc  convergent  nozzle  installed. 
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of  Characterisitcs  Program  based  on  the  appropriate  NPR.  No 
inviscid/inviacid  iteration  was  performed,  that  is  to  say, 
the  external  pressure  was  assumed  to  be  a constant  ambient 
value  with  no  longitudinal  variation.  Thus,  it  must  be 
realized  that  the  plume  shapes  used  in  the  analysis  were  most 
likely  too  large  by  a small  amount,  as  the  increased  pressure 
in  the  compression  at  the  end  of  the  boattail,  relative  to 
tunnel  ambient  pressure,  would  tend  to  decrease  the  plume  size 
in  that  area.  Recovery  to  near  ambient  tunnel  conditions  was, 
however,  fairly  rapid. 

To  provide  the  velocity  ratio,  for  use  in  the 

mixing  program,  DELPL,  the  velocity  of  the  exhaust  flow  under- 
going isentropic  expansion  to  tunnel  ambient  pressure  was  used, 
along  with  the  velocity  of  the  wind  tunnel  evaluated  also 
under  ambient  conditions.  To  furnish  the  appropriate  con- 
ditions at  the  end  of  the  nozzle,  including  the  internal 
boundary  displacement  thickness  at  this  point,  a one-dimen- 
sional isentropic  analysis  was  used  based  on  plenum  and  nozzle 
geometries  and  known  plenum  conditions.  The  resulting  boundary 
layer  growth,  calculated  using  the  Bartz  program,  is  shown  in 
Figure  26.  Also  shown  in  this  figure  for  comparison  purposes 
is  the  result  of  the  same  calculation  done  with  the  recently 
available  axisymmetric  version  of  the  Kuhn-Nielsen  Separated 
Turbulent  Boundary  Layer  Program  (77)  , utilizing  its  boundary 
layer  option.  The  final  value  of  the  displacement  thickness 
of  the  nozzle  was  the  value  used  as  input  to  the  mixing  pro- 
gram. The  corresponding  value  of  the  external  boundary  layer 
displacement  thickness  at  this  point  was  provided  by  the  most 
recent  calculation  of  the  viscous  portion  of  the  CVI  Program, 
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Having  specified  all  of  the  requirements  of  the 
mixing  program,  a number  of  calpulations  were  made  with  the 
CVI  Program,  including  DELPL,  for  different  values  of  the 
mixing  parameter,  a.  As  shown  in  Figures  27a  and  b,  for 
NPR's  of  two  and  four,  respectively,  the  circular  arc  boat- 
tail  was  analyzed  with  the  appropriate  plume  shape  for  the 
"solid-body  plume"  case  and  for  the  mixing  cases  with  the 
following  values  of  o,  which  were  described  in  Chapter  III, 
page  79; 

(a) 

(b) 

(c)  cr^/10 

It  is  emphasized  that  the  plume  shape  was  held 
constant  during  this  analysis  for  purposes  of  isolating  the 
entrainment  effect.  The  value  of  u may,  in  actual  cases, 
never  be  as  small  as  case  (c) , o^/lO,  since  part  of  the  6* 
correction  would  be  taken  up  by  the  actual  changing  plume 
shape.  The  choice  of  c^/lO  was  arbitrary  and  was  used  to 
note  the  sensitivity  of  the  solution  to  changes  in  c. 

For  the  two-stream  value  of  c?  in  case  (a)  , the 

difference  between  the  C distribution  and  that  of  the  solid 

P 

body  plume  was  so  small  as  to  be  negligible.  There  is  a 
slight  improvement  when  the  one-stream  value  of  a , case  (b) , 
is  used  particularly  in  the  NPR  = 4 case,  where  the  CVI  Pro- 
gram prediction  agrees  well  with  the  data  except  at  the  last 
two  data  points,  which  may  be  an  indication  of  localized 
separation  in  that  area.  A small  separated  area  was  also  a 
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Figure  27.  Longitudinal  pressure  coefficient  distribution  for  the 
NASA  Langley  2./d  = 1.768  circular  arc  boattail  with  a real 


Figure  27.  (Continued) 
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chcirs-Ctef  is  tic  of  th©  d3.tB.  in  the  solid  sting  case. 

It  was  arbitrarily  decided  to  run  these  cases  with 
an  extremely  low  value  of  c , case  (c) , which  would  be  indic- 
ative of  very  high  mixing.  As  can  be  seen  from  Figure  27, 
this  improves  the  results  considerably  for  the  KPR  = 2 case, 
yet  causes  a large  over-correction  in  the  NPR  = 4 case.  This 
is  a result  of  the  much  higher  entrainment  rate  at  a given  a 
value  for  the  NPR  = 4 case,  which  has  a value  of  0.585, 
compared  to  the  rate  for  the  NPR  = 2 case,  which  has  a <t>j^ 
value  of  0.790,  due  to  the  difference  in  The  large  en- 

trainment rate  produces  a much  larger  6*  in  the  NPR  = 4 case. 

It  may  be  seen  from  these  cases  that  there  is 

great  room  for  improvement  in  the  choice  of  a for  each  point 

in  the  mixing  calculation  for  a given  configuration.  But, 

qualitatively,  it  has  also  been  shown  herein  that  applying 

the  6*  correction  at  points  downstream  of  the  nozzle  exit  is 
e 

sufficient  to  influence  correctly  the  pressures  over  the 
entire  afterbody  and,  thus,  that  the  displacement  thickness 
correction  method  is  a valid  approach  to  the  plume  entrainment 
problem. 

An  example  of  the  typical  behavior  of  the  6*, 
and  y in  the  area  of  the  boattail  and  plume  is  shown  in 
Figure  28. 

AGARD  15-deq  boattail . Another  axi symmetric  "cold" 
flow  exhaust  plume  model  was  also  analyzed  with  the  CVI  Pro 
gram  and  compared  with  data.  This  model  was  the  AGARD  15-deg 
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Boattail  (23) , and  it  was  chosen  because  of  its  geometry  and 
the  NPR's  at  which  it  was  tested.  As  can  be  seen  from 
Figure  29,  the  boattail  is  not  only  steeper  near  the  exhaust 
plane  than  the  Langley  = 1.768  circular  arc  boattail, 

but  it  has  a slightly  lower  closure  ratio  ‘ addi- 

tion, the  NPR's  at  which  the  AGARD  15-deg  boattail  was  tested 
covered  a greater  range  than  those  for  the  Langley  boattail, 
three  to  seven  in  this  case  versus  two  to  four  in  the  former 
case.  It  was  felt  that  the  combination  of  these  two  effects 
would  provide  a more  challenging  case  than  the  previous  one. 

The  CVI  Program  Analysis  was  carried  out  at  the  two 
NPR's  at  a Mach  number  of  0.8.  The  results  were  compared 
with  data  from  Reference  (23)  , in  which  the  model  was  tested 
in  the  AEDC  16-foot  Transonic  wind  Tunnel.  The  internal 
plumbing  arrangements  to  provide  high  pressure  air  to  the  jet 
exhaust  were  very  similar  to  those  of  the  Langley  tests,  thus 
an  isentropic  one— dimensional  analysis  was  utilized  to  obtain 
conditions  in  the  plenum  and  at  the  nozzle  exit.  The  methods 
of  estimating  the  necessary  boundary  layer  and  velocity 
parameters  for  use  in  the  mixing  program  were  also  used 
exactly  as  in  the  previous  subchapter. 

The  comparisons  between  the  theory  and  the  data  are 
shown  in  Figures  30a  and  b.  General  agreement  between  the  two 
is  fair  in  both  cases,  but  some  points  must  be  made  concerning 
the  comparisons.  Firstly,  the  grid  points  at  the  shoulder  of 
the  boattail,  where  the  flow  is  expanding  the  most,  miss  the 
point  of  minimum  pressure.  This  particular  misplacement  of 
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Figure  29,  Comparison  of  boattail  and  plume  simulator  shapes 
between  the  NASA  Langley  = 1.768  circular  arc  boattail 

and  the  AGARD  15-deg  boattail. 
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a.  NPR  = 3 

Figure  30.  Longitudinal  pressure  coefficient  distribution  for  the 
AGARD  15-deg  boattail. 
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points  could  not  have  been  predicted,  nor  could  it  have  been 
recognized  without  the  comparison  with  the  data  shown  in  this 
figure.  It  would  be  suggested,  then,  that  future  analyses 
cluster  the  points  in  this  area,  and  that  once  a successful 
analysis  has  been  completed,  the  grid  points  be  shifted 
slightly  to  check  for  sharp  gradients  near  the  shoulder. 
Secondly,  there  appears  to  be  a rather  large  separated  region 
near  the  end  of  the  boattail.  This  is  especially  noticeable 
in  the  NPR  = 7 case.  This  invalidates  the  results  of  the  CVI 
Program  in  this  area  of  the  analysis. 

Thirdly,  the  results  of  changing  the  value  of  a from 
CTjj  to  oj-  and  then  to  aj./10  are  more  consistent  in  these  cases 
than  in  the  Langley  circular  arc  boattail  case-  This  may  be 
due  to  the  fact  that  the  values  in  this  case  were  closer 
than  in  the  previous  case,  0.653  and  0.530  versus  0.790  and 
0.585.  Thus,  the  difference  in  entrainment  characteristics 
and  hence  6*  variations  were  not  as  pronounced  in  this  case 

as  in  the  previous  case. 

The  analysis  of  the  AGARD  15-deg  boattail  thus 
points  out  the  qualitative  ability  of  the  CVI  Program  and  the 
mixing  model  to  predict  real  plume  flows,  but  it  is  apparent 
that  further  work  must  be  done  on  the  determination  of  the 
longitudinal  variation  of  o at  each  iteration,  and  that  the 
important  and  complex  interaction  between  the  plume,  external 
inviscid  flow,  and  separated  boundary  layer  must  be  studied 
in  detail. 
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CHAPTER  V 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  objective  of  this  study  was  to  develop  a method 
for  prediction  of  the  transonic  flow  about  an  arbitrary  axi- 
symmetric  body.  The  impetus  for  the  effort  was  the  desire  to 
compute  the  flow  over  nacelle  afterbodies  for  purposes  of 
drag  estimation,  but  this  objective  was  quickly  extended  to 
include  other  axi symmetric  bodies,  such  as  equivalent  bodies 
of  revolution.  Such  flows  may  exhibit  strong  interactions 
between  the  external  inviscid  flow  and  the  viscous  flow  near 
the  body,  which  may  include  large  regions  of  locally  sonic 
flow,  near-normal  shock  waves,  and  extensive  separated  areas. 
If  the  axisymmetric  body  also  includes  a jet  exhaust  plume, 
the  interactions  are  even  more  complex,  including  now  the 
inviscid  flow  of  the  plume  and  the  turbulent  mixing  between 
the  plume  flow  and  the  external  flow. 

Many  analytic  methods  for  solving  problems  of  this 
sort  were  examined.  Following  this  effort,  a combined 
viscous/inviscid  iterative  digital  computer  program  was 
developed  which  served  as  a means  for  obtaining  the  flow  pre- 
dictions, including  the  effect  of  boundary  layer,  for  com- 
parison with  various  data.  Finally,  a turbulent  mixing 
analysis,  based  on  the  Chapman-Korst  mixing  theory,  was  used 
to  determine  the  effect  of  plume  entrainment  on  the  afterbody 
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pressures.  The  conclusions  and  recommendations  below  are  the 
result  of  all  of  the  efforts  outlined  above. 

1.  Examining  contemporary  works  revealed  that  the 
viscous/inviscid  approach  to  problems  of  this 
sort  was  by  far  the  most  popular.  One  interest- 
ing approach  used  the  Navier-Stokes  equations 

in  the  af terbody/sting  area. 

2.  The  comparison  of  the  results  of  the  method  of 
analysis  developed  herein  to  data  leads  to 
certain  conclusions  with  respect  to  the  physical 
phenomena  involved; 

(a)  The  predictions  agree  well  with  the  data 
if  the  flow  is  unseparated.  This  re- 
stricts the  angle  of  the  end  of  the  boat- 
tail  to  approximately  20  degrees , with 
larger  angles  tending  to  cause  separation 
of  the  flow  at  almost  all  Mach  numbers. 

(b)  Care  must  be  taken  in  analyzing  the 
results  of  the  method  at  the  upper  tran- 
sonic Mach  numbers,  where  the  separation 
problem  is  greatly  intensified. 

(c)  The  slight  overcompression  predicted  at 
the  junction  of  the  afterbody  and  solid 
sting  is  there  even  in  otherwise  well 
behaved  cases.  This  may  indicate  a small 
separated  area  or  an  inaccuracy  in  the  value 
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of  (5*  at  this  point  for  all  configurations 
and  Mach  numbers,  even  quite  moderate  ones. 
Thus,  slight  inaccuracies  at  this  point  are 
to  be  expected  even  if  the  extreme  condition 
of  (a)  and  (b)  are  avoided. 

(d)  The  standing  shock  developing  at  the 
afterbody  shoulder  which  increases  in 
strength  and  moves  backward  with  increas- 
ing Mach  number,  is  predicted  quite 
accurately  by  this  method. 

(e)  For  most  bodies,  the  boundary  layer  effect 
is  not  important  upstream  of  the  boattail. 
This  will  not  be  true  for  bodies  with 
exotic  forebody  curvatures,  however. 

Concerning  the  CVI  Program: 

(a)  Relaxation  of  the  addition  of  the  boundary 
layer  displacement  thickness  may  increase 
the  rate  of  convergence  or,  in  some  cases, 
make  the  difference  between  convergence  or 
divergence . 

(b)  The  South-Jeiraeson  Transonic  Program 
(RAXBOD)  is  well  suited  for  use  in  the  CVI 
Program.  Care  should  be  taken,  however, 
in  the  results  of  some  of  the  calculations 
for  bodies  with  high  body  radius  gradients 
at  high  subsonic  Mach  numbers,  particularly 
if  the  results  vary  noticeably  in  form  from 
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from  results  at  lower  Mach  numbers  (see 
Appendix  C) . 

(c)  The  Bartz  boundary  layer  program  has  the 
capability  of  continuing  the  solution 
through  areas  of  high  pressure  gradients, 
which  makes  it  valuable  in  the  initial 
stages  of  the  viscous/inviscid  iteration. 
This  same  characteristic  will  sometimes 
produce  converged  solutions  in  the  CVI 
Program  even  if  the  flow  is  actually 
separated.  Any  unusual  behavior  of  the 
Bartz  program,  particularly  any  unusual 
variation  of  the  displacement  thickness 
calculation,  must  be  considered  suspect. 

It  is  recommended  that  this  program  be 
replaced  by  a more  accurate  boundary  layer 
program  in  the  latter  stages  of  the  CVI 
Program  iterations, 

(d)  It  is  recommended  that  a method  of  char- 
acteristics program  be  incorporated  into 
the  CVI  Program  to  enable  calculation  of 
a real  plume.  This  would  provide  an  auto- 
mated inviscid/inviscid  iteration  capa- 
bility in  addition  to  the  viscous/inviscid 
iteration  capability  presently  existing. 
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4.  The  application  of  the  Chapman-Korst  mixing 
theory  to  the  modification  of  the  plume  bound- 
ary has  the  advantage  of  easy  implementation  as 
well  as  being  physically  satisfying.  The 
derivation  of  the  equations,  although  restricted 
to  two-dimensional  isobaric  flow,  included 
terms  accounting  for  the  variation  from  the 
fully  developed  similar  solution  of  Korst  and 
Chow  (51)  due  to  initial  boundary  layer  shape. 

The  displacement  thickness  correction,  5*,  due 
to  plume  entrainment  was  applied  directly  to  the 
boundary  layer  which  was  calculated  assuming  the 
plume  to  be,  at  first,  a solid  body.  Certain 
modifications  of  the  theory,  however,  are 
recommended:  (a)  the  basic  equations  should  be 

rederived  for  the  axisymmetric  case  and  a longi- 
tudinal pressure  distribution,  and  (b)  the  mix- 
ing parameter,  c , should  be  calculated  on  the 
basis  of  conservation  of  momentum  considerations 
at  each  longitudinal  point,  as  in  the  work  of 
Bauer  and  Matz  (47) , rather  than  computed  on 

the  basis  of  previous  work,  as  is  currently  done. 
The  values  of  o should  be  recalculated  at  least 
everytime  that  the  plume  shape  is  redefined. 

5.  It  is  felt  that  the  highly  complex  problem  of 
turbulent  separated  flow  with  plume  entrainment 
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effects  is  of  sufficient  practical  importance 
to  be  investigated  in  future  efforts.  In 
particular,  it  is  recommended  that  the  recently 
developed  Kuhn  and  Nielsen  Axisymmetric 
Turbulent  Separated  Boundary  Layer  Program  (77) 
be  modified  by  the  addition  of  a plume  entrain- 
ment segment  of  the  program  based  on  the 
Chapman-Korst  mixing  theory  similar  to  the  work 
contained  herein. 
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APPENDIX  A 


DEVELOPMENT  OF  AN  APPROXIMATE  CLOSED  FORM  SOLUTION 
FOR  THE  TRANSIENT  INITIAL  BOUNDARY  LAYER  INTEGRALS 


A method  of  approximating  the  integrals  encountered 
in  the  expression  for  <)),  Equation  (8),  page  60,  was  suggested 
to  the  author  by  Whitfield  <78}  and  is  based  on  the  series 
expansion  work  of  Chapman  and  Korst  (66)  . Implementation  of 
this  series  approximation  has  decreased  the  CPU  time  involved 
in  the  evaluation  of  these  integrals  using  Gaussian  quadrature 


by  a factor  of  five  with  negligible  loss  in  accuracy,  even  at 

4 

extremely  high  values  of  on  the  order  of  10  . 

Basically,  the  method  involves  the  realization  that 


evaluation  of  the  integral 


(A-1) 


depends  strongly  on  the  extremely  peaked  nature  of  the  ex- 
potential function.  As  such,  the  boundary  layer  profile  in 
the  integrand  may  be  expanded  in  a Taylor  series  about  the 
maximum  value  of  the  integrand,  which  occurs  at 

^max 

The  positive  root  of  the  differentiation  has  been 
rejected  since  it  always  falls  above  the  upper  limits  of 
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integration.  For  all  intergration  limit  values  of  n and  rip» 

expansion  about  this  maxiirium  is  valid,  with  one  exception. 

Only  if  the  lower  limit  of  integration  exceeds  the  value  of 

8 is  this  not  possible,  in  which  case  the  integrand  may  be 
max 

expanded  about  the  lower  limit,  n “T1„* 

Given  I-  as  defined  in  Equation  (A-l) , the  velocity 

& 3 

profile  of  the  boundary  layer  may  be  expanded  in  a Taylor 
series  for  the  case  n - Hp  < ^max'  follows; 


gtsl  = 


= g ( 3 ) + 


dq  (B  ) 6-B 

SB “TT 


thus 


■2a 


/ 1-^ 
n-n„  \ 


, 6-B 

1 max 

m n-6 


max 


e-e  dB 


or 


/T~  m(n-B„.„) 

max 


[m.  (l-m)e^ax]  / 


n -n. 


e-e  dB 


/ 


h-n. 


0 6®  dB 


which  becomes,  after  integration. 
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"2a  - 


n-e  ,1/in  / 

max  ' ' 


2/?m(n-6jj^) 


[mn+(l-m)e  ] /it  erf  n - erf(ri-n_) 

In&X  I P 


2 , ,2 


If  n - n ^ , then  the  corresponding  expression  is; 

p IllclX 


■2a 


2/tt  mn 


<[n  +(m-l)Ti  ]/tt 


erf  n - erf  (n-np)J 


Similarly,  if  the  upper  initial  boundary  layer  contribution 
is  considered,  that  is  to  say,  the  integral  involving  '^2h' 
then  the  approximate  closed  form  solution  may  be  written  as 


■2b 


& n-3 

a max 

^ n 


ll/m 


2/i 


k{ 


[mn+  (1-m) 


e£f  |n  + 3^ 


- erf  nj  + e 


6. 

- 2 
n ' 


- e 


for  the  case  where  n + n > P , and 

5 p max' 
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■2b 


Z/tt  m -j—  n 


T-  "p 
a 


r 1 

n-(m-l)  3-  n 

L a 

/? 

erf 

n+  y-  n 


- erf  n 


“ ^ 2 
+ e a ^ 


b 

for  the  case  where  n + ^ n - 

a ^ 
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APPENDIX  B 

LISTING  OF  THE  PLUME  ENTRAINMENT  PROGRAM,  DELPL 

The  DELPL  Program,  as  described  in  Chapter  III,  and 
its  associated  subroutines  are  inserted  directly  into  the  CVI 
Program  when  an  afterbody  with  a real  plume  is  analyzed. 

DELPL  obtains  its  required  input  parameters,  which  are  listed 
in  Table  VI,  page  68,  from  a calling  subroutine  in  the  CVI 
Program  through  its  argument  list  and  returns  the  values  of 
the  displacement  thickness  corrections  due  to  plume  entrain- 
ment in  a like  manner.  The  numerical  integrations  are  per- 
formed by  the  subroutine  GSS96,  which  is  called  by  DELPL  at 
several  locations.  GSS96  in  turn  calls  the  function  F, 
which  evaluates  the  velocity  function,  , and  provides  the 
numerical  values  of  the  integrands  of  the  integrals 

I.,(n  ),  and  other  integrals  used  in  the  calculation  of  the 

2 Ra 

stagnation  temperature  correction. 

The  integrals  appearing  in  the  basic  relationship 
for  4),  Equation  (8),  page  60,  are  evaluated  in  function  F 
using  the  approximations  described  in  Appendix  A at  a large 
saving  in  CPU  time. 

The  DELPL  program  is  relatively  fast  compared  to 
the  CVI  Program,  with  twenty-five  complete  6*  calculations, 
a typical  number,  requiring  about  ten  CPU  seconds  on  the 
IBM  370/165.  For  a typical  viscous/inviscid  solution  of 
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twenty  complete  iterations,  the  increase  in  CPU  time  due  to 
the  entrainment  calculations  is  approximately  10%. 
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SUenOtlTiNL  DtLPL  {X>fQ«TUAtTUH«C:P*«LtA<DSTHAtUSTPUiE.NiMi  IWH.OSTPEl 
IMPLICIT  SEAL»BlA-M,0-Z) 

QIMtMSlON  X 1100)  lOSTREdOS)  f JMN(3} 

OIHCNSION  ETAUUMIlOll  fFirUHdOLl 

COMMON  /PAPAM/£TAtaETAiSQP1N<FFtt«ETAP400A.01JR<TT»CAZ>P<iP 

COMMON/ XDUMH/XOOM 

COMMON/TSeBC/F5(JliC.£TAL«ETA« 

SUPIN».Si6AlB<»&B3bDO 

FSLeC>>.87A6?ST3D0 

» MECEJXED  IN  IhCHESi  THEN  CONVERTED  TO  FEET  FO»  CALCULATIONS. 
«•»  RETUHNtO  IN  INCHES. 

DO  lOA  l*liN 

IFIX  dJ  .LE.O.UO)  Xd)=I.D-3 
Xd)=Xd)/12.D0 
4 CONTINUE 

COMPUTE  CROCCO  NO.  SOUAKED 

CAZ»ua»«2/  [?,U0*CPA»T(JA) 

***  COMPUTE  OTHER  NEEDED  PARAMS, 

DA= I I .UO*EN) •DSTRA 

Dd=l l.DO«EN)*DSTRS 

Tma*EN/(2.O0*EN) *OSTRA 

THB»EN/)2,D0*ENI *OSTRd 

DDAsUA 

DOBsDB 

TT=T0B/T0a 

PWP>I .DO/EN 

FFb=Fa 

Fa2=FH»»2 

RRsd .D0-CA2)/<Tr-CA2«F«2) 

C 

C INITIAL  OUTPUT 

C 

MRITElfti 110) 

110  format (•0'.Tfl*‘CA2».T23i'OA'»T38.iOd'.T&3.»THA'.T&8.'THH*.r«3. 

« • r T • . T 9d  . ‘PXIR  • . T 1 1 3 . • HR  ' ) 

MHIT£C6.1JS)CA2.OA.0e.THA.TMBi TTiPdR.RR 
IDS  FOHMATIIH  40  dPtlS.a)  f/l 

IF  duR  d ) .EO.  1)  ■fllTEIb.lOl) 

101  FORMATdl'.TbT.'COMPUTE  0»£*4//) 

IF  [iNRd)  .EO.l)  mRITE  (E.ioai 

10  2 FORMAT  1 TR.<XII)<4T24  4 'SIOMA  tiTAO.’ETAP'.TST.'Il'tTTS.'l^'.TAR. 

**IElJ*.T105.*ETAM«,T12tl.*DSTHE«/) 

OO  ICO  i>l.N 
C 

C •••  N IS  NOMHER  OF  X'S 
C 

xOuM«12,DO»x ( I) 

C 

C •«*  COMPUTE  SIGMA.ETAP.  ETRA,  AND  ETRd 
C 

SIG  >12.DO*2.TSeOO*D5URT )S.U0>CA2/ ( 1 .D0-CA2) ) 

SIB»d.U0*FB) /d,00-Fd)*SIG 


US  THE 


C 
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51G0*SI6 

SIG«SIG« (X ( 1 1 / /3b.q21U0P •• 

IF (SIG.GT.SIGO)  SIG=SI&0 
ETAP  s SIG<DA/X  1 1) 

ETBAi»ETAP*a.OO 

ETR»"-Gb«E7AP/DA-A.O0 

C 

c »|>»  iteration  for  stag.  temp,  correction  integral  limits 

c 

ETAH=ETHA 
OETAHotTRA/a.Dl 
ETAL»eiHB 
UtTAL=-tTRe/2.01 
231  ETAM»ETAH-OeTAH 
FETAHbF (ETAH,3» 

IF  1FETAM.LT-9.99'#D-1)  GO  TO  230 
GO  TO  231 

230  IF (DETAh.LT.1 .0-b“£TAPl  GO  TO  232 
ETAH*ETAH*oeTAH 
DETAMaOtTAH/1 0 .DO 


GO  TO  231 

232  ETAL»ETAL*DETAL 
FETAL^F (ETAL»3t 

IFiFETAL.GT. 1 .OOOlOOoFa.OR.F ETAL.LT.S.V^RO-l^FBl 

GO  TO  232 

233  IF [OETAL.LT.1.0-5*DB/UA»ETAPI  GO  TO  23* 
£TAL=ETAL-DETAL 

OET  AL"17ETAL/10.U0 
GO  TO  232 
23*  CONTINUE 
C 

c COMPUTE  Integrals  for  stag,  temp,  corhection 
c 

call  GSS96(A*ETAL.£TAH,PN1ERH 
CALL  GSS96l3»ETALiETAM,PMlACt) 

IF ll.EO.n  FSCREF-tPKlEHF-PHlACT) /PHitRF 
FSBCI® (PHlERF-PNlACTl/PrtlERF 
F5uaC»F&eCl/FSC»tF 


GO  TO  233 


C 

c COMPUTE  IHETARA)  ANO  I2ltIARAl 

C 

SUMi laO .DO 

SUNI2>0.D0 

NLN*1 

XNLMahLN 

00  220  NLal.NLM 

XNLBNL 

ErRL»ETRB» IETHA-ETR91 /XNLR» t XNL- 1 . OO > 

ETRMaETRB*  tE TRA-ETRB > / XNLN*XNL 
CALL  GSSRGd.tTRL.ETHH.Sll) 
call  GSS96(2*ETRL«eTfiH.5I2> 

SUM] l-SUMIt^Sl 1 
SUMI2aSUM12*SI2 
220  CONTINUE 

SIP  = tfRB«(l.D0-CA2l«FiJ/ITT-CA2»Fd2l 

SIl»SUHll.SIP 

SI2»SUMl2»SIP*Fb 

SIlCHK»tTRA-£TAP/OA® (USTRA»HR»Fd*D5IRai 
SI  2CMXatTRA»ET  AP/DA®  (DSTRA*TMA»HR*F  B2®  (USTRb*  Trta  I I 
210  FORMATdH  ,SdPE15,8l) 

C 

c ••O  COMPUTE  II(ETAJ) 

C 

SllEJ»(SI  l-S12-tTAP»<fiR“FH2‘’TH6*THAI  /DAI  / d .OO-Fd! 

C 

c ®««  compute  ETAM 
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ETAH  =»  LIRA  - ( S I 2-FB*S  1 I ) / U . UO-F  6) 

ETAM  ■ ETAM  - eTAP*(TMU»RR«FH2  ♦ ThA  ♦ J5TRA»(1.D0-  FBI)/ 

» (0A»(1 .DO-  FH) > 

»•*  COMPUTE  DSTREin 

DSTnE(  II  = 0A*(ETAM.SI  IE  J«  I T T-CA2»Fb2  ) / 1 F (1.0n-CA2)n  /FI  AP-OSTHB 

OSFRE  CONVERTED  10  INCHES.  X CONVERTED  BACK  TO  iNCntb. 

OSTRE (1) ■12.0  0*0 S THE  111 
X(II >12.D0>X  in 

IFIluHIlI.EO.n  iiRITE(fc«lQl)xiI)»SI6i£TAP,SlUS12iSIlEJ.tlA*‘. 


•OSTfiE(l) 

103  F0RMAT<6(1PE16.6J  I 

100  CONTINUE 
RETURN 
END 

SUBROUTINE  GSS96  IK. A. B. El  AG  ] 

IMPLICIT  REAL«B(A-H,0-7)  Ab  ^ 

COMMON  /CARE/  NEAT  Ab  3 

DIMENSION  X(*a).  mlABI  AC  A 

data  N/*R/  AG  5 


DATA  X/.DlfePTbTAARAREOP^S.  .OA8ei2'»a!5l360AR70.  .OH1?97*R5*6A4?S*0.  .1  AG  f, 
UaGRSMbOl  J 066S7,  .14597311  *6Sfc896T».  I 7809ft8923N7bI8S..2l  0 031  31fl4N0S  AG  7 
2871 • .2* 174 3 1581 638399.. 273 19881 2591 0489.. 30438494435*4982. .3382085  AG  8 

3228928283. . 385898881 4 7231 35.. 395 7976498289084. .425*789884073005. .4  Ab  9 

4547094221877428. . 4834579739205962..5 11894 17 7 154887 7.. 539388108 3243  AG  10 

5572. . 5885104105813971. .59303236*7775719. .8189258401254885. .8441634  AG  11 

8037849870..  668718310 0439160. .6925645368421715.. 7158768123*89874. • 7 Ab  12 
7380  30  6437*4  40  0 0, .759  8023*1  1766*  7 2. .7803690*386  7*3  31  ,.8003087**1391  AG  13 
8*0  7,,  81 9*00  31  073793 15..  837823511 228 1 869..a54959D3  3 *3*6(1  13.,a  7 1 3885  1* 

9059092963.. 8B6e9*51 7*02*203. .901*606353158521 ..9 1507 1*2312 08 981. .9  AG  15 
»27712*567223085,.9393T0339T527551,.950032 7 I 778**37*. .959888291**87  AG  18 
**23,.98B328a28*e328*0, , 9759J91 7*585 1 38* .. 9825 1 72 63563 0 1 *6 . .98805* I AG  IT 

12632962 36..  9925* 39 0 03 2376 24, , 9 9598 1 8*29872 09 l ,. 998 364 3 75  883  I 8 1 5 . .9  AG  10 
1996895038832307/  Ad  19 

data  m/. 032550614*9236318. .03251811871386891., 0324471637140642*. .0  AG  20 
i 323*38  225885759  I,. 0 3220620* 79*03022,. 0 32034*5623 19926*,. 031 028758  9 *5  21 
294*1 099.. 031 589330770727  I*. .031316*256968813*. .03101033250631383,.  AG  22 

303067137612366912. .  0302999 1 5*2082 75 7 02989634* 1 3632630. . 029*6 1 089  A3  23 
*95816790,  .02899*61  *15055522..02R*97*U  06508536.  .0279  7000  76168*831  . AG  ?* 
5. 027* I 2 9627 260 2922,. 02682686672559 17*.. 0262123*0 7 35672*0,. 02557003  AG  25 
68005349  35.  . 02*90  0833222*0359  0 2*20*8*  I 792  36486 ..  0 23*0  3 3990  ati92820  AU  26 

7..  02273706985832935.. 02 19666***38 7**32. .021 17293 9092 191 29., 0203567  AG  27 
89715*33329, .019519001 1*01*500,.018660679627*H*5, .0177025023160*52  AG  28 

93. . 0 1688547966*2451 b..0l 59705629025622 7, ,015038 721 02699*93,.01 *090  Ad  29 
«9*177231*85..01312S22956696 156.. 01 215160*67108029,. 01 1182102099036  Ab  30 
**e. .010160  770535008*1  I .0091*8671230  783383. ,0081266  78925898756. .007  Ab  31 
*096*70  791  I 53859. ,0  080585*550*2  35957 .,0050 1*20 27*292 751 *,.00398 *55*  AG  .32 
*338***883, .00291073181793*9*4, .001 B5398OT8094592O..OOO798792O65552  AG  33 


*0122/  Ad  34 

e=o,ooo  »e  35 

HR=0.5*IB-A!  AG  36 

FMIO*a,5*(A*8)  AG  37 

DO  1 .J»1,N  AG  38 

«TsF«ID.HR»X I J)  AG  39 

XS»FMI0-HH»X IJI  AG  *0 

e«E**  ( J)  *IF  (XS.M  .F  (XT  ,K)  I AG  *2 

1 CONTINUE  AG  43 

E»MR»E  *U  *4 

RETURN  Ab  *5 

EnU  Ab  *6- 

F UNCT ION  F I X. I I 
implicit  REAL'S  I4-H,Q-Zl 
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COMMON  tTA»BEIA«SUPlNtPB«ETAP.UAiUB>TrtC:A2tPKH 

COMMON  i'XOJKM/XDUm 

COMMON/E  Suec/F  5UHC t ETAL  »ET  AM 
£TA»* 


*»»  51VEN  ElAt  FIND  HMl 


PH1=0.SU0*I1 •U0*FB) 

PMZ»0.5l)0*D£HF(tIA-ETAPI 

PH3»-0  «bOO*Frt*OtHF [E T A*oe*E T AP/OA ) 

«•«  CALCULATION  OF  INTEGRALS  IN  VELOCITY  FUNCTION  BY  SEHILS  EXPANSION 


3 


A 


b 


6 


TEPHAbET A>ETAP 

TERMb»ETA*DB/l>A*tTAP 

DLA30.U0 

ULAsl.UO 

DL5=-0U/DA 

UL5a0>O0 

IF<ETA.GT.-A.OOiANO.ETA.LT.ETAP»A,OOt  GO  TO  3 

PHAiO.OO 

eO  TO  A 

CONTINUE 

•BG1»ETA»*2 

ARGE»(ttA-£TAP)**2 

IF  IAft61.6T.Ub.OOI  AftGl«Jl5.00 

IF IA»62.6T.llb.D0l  AfteZ»llb.O0 

9ETAMA«tTA/2.l)0-05QBT  ( ( ETA/2  .DO ) .00  I 

TCa-DEBF (ETAI-OEHF(TERMA) 

TOA=DEXP  < -ARG 1 ) *OEXP (-AH62 I 
IF < TEHHA.OE.BETAMAI  60  TO  20 

TAA«<  (ETA-HETAMA)/ETAP)«»Pli(«/2.O0»SOPI«-**NP/I£TA-eETAM4) 

TBAsIETA/RKR*  U .OO-I.OO/PPR)  *bETAMAI  /SCjPlN 
GO  TO  21 
20  CONTINUE 

TAA«l.DO/2.DO»SOP!N«PPB/eTAP 

TBA»(EIA*U.OO/P»B-I.OO)»ETAP)/SQP1N 


21  continue 

PMA«TAA*tTBA*tCA*TOA) 

continue 

lF{£TA.LT.A.OO.AND.eTA.6T.»OB/OA»ETAP-A.UOt  GO  TO  5 

PHbaO.OO 

GO  TO  b 

CONTINUE 

APQla IETA«Pa/UA*ETAP)**2 
Afta2aETA»A2 

lF(APGl,GT,Ub,D01  AP61«U5.l)0 
lF(APG2.6T.Ub.OOI  APG2-U6.U0 

6ETAM5»tTA/2.00«OSOPT( IETA/2.UO)**2*PWB/2.O0l 
TC5=0EHF (TERMbI -OERF (ETA» 

TObaOEXP I >ARG 1 1 -DEXP  <-AM62 1 
lF<TERMB,LE.8ETAH5»  60  TO  30 

TAS=  t«OA/OB*(ETA-0ETAMbl/ETAPI  a*pbp/2.D0»SQP1N*®pR'' (E  TA-dET  AMSi  • 
rab»(£TA/PvB» ( 1 ,00-1 .OO/PdP) •etTAHbl/SOHIN 
GO  TO  31 
3D  CONTINUE 

T A5  = -FH/2,00*S(lPlN»PdP/tTAP/DB*UA 

T8b=  (ETA-Cl.OO/PpP-l.OO)  ‘OB/IJAaETAPJ  /SOP IN 


F« 


31  CONTINUE 

PM6«TA5A(TB5«TCb«T05) 

CONTINUE 

PMI=PMl*PM2*PM3*PHA»PMb 


C 

C ADJUSTMENT  TO  STAG,  TEMP.  PAPAMETEB.  LAHOA 


C 


B]6LAM«TT* ( 1 .UO-BMI 1/(1 ,O0-F  B) • (PMI-F0) / ( 1 .UO-FB  > 
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IMET*,(»C.ET»t'l  faO  TO  900 

IF  <ETA.U€,-OH/OA«ET4f>)  SO  TO  900 

0=7.390-2 

c=c«fsubc 

B[6L4W  = BIGLAM-C 

IF (ETA. G€, 0.00)  SIGLAM  = BI 6l4«*C» ETA/t T AP 
IF  (ETA.UE,  0.001  UlGLAK=BieLA«-C«tTA/tI4P*i>A/0H 
900  CONTINUE 

FUNC«I 1 .O0-CAe)/(BlGLAM-C42»PHI»*2) *PHI 

If IETA.lT. 0.001  SO  TO  10 
P<i]aL=l.O0 

If IETA.LT. ETAP)  PHl9L»(tTA/tTAPI«*PWR 
SO  TO  IJ 
10  PHiaLsFB 

IF (ETA.GT.-DB/UAeETAP)  PM  I BL  = FBe ( -OA/DUeE I A/ETAP ) «ePBR 
1]  CONTINUE 

GO  TO  II  iZ.AO.Al  )'i  1 

*•*  F FOP  1 1 (ETAPAI 

] CONTINUE 
FaFUNC 

100  FONHATIlM  , lA.  10 < IPElZ.bl 1 
KETURN 

• ■»o  F FOR  IZtETARAI 

2 CONTINUE 
F=FUNC«PHI 

return 

•••  F IS  VELOCITY  PATIO  FUNCTION  FOR  ACTUAL  FlOB 

40  CONTINUE 
F»PHI 
RETURN 

•••  F 15  VELOCITY  RATIO  FUNCTION  FOR  FULLY  UEvElORED  F lO« 

41  CONTINUE 

I F lETA.sE.O.DO)  ETAT«P=1 .aSOOeETA/E T AM 
IF  IETA.LT, 0.00)  ETATHP= 1 .B3U0e£1 A/ETAL 
F = O.SDOe I I ,DOeFb*DERF IE TATMP) -FB»0ERF (ETATMP)  ) 

RETURN 

end 
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APPENDIX  C 

ANALYSIS  OF  THE  CBA  BOATTAIL  - AN  ANOMALOUS  CASE 

An  analysis,  using  the  CVI  Program,  of  a boattail 
similarly  analyzed  by  Chow,  Bober,  and  Anderson  {9) , resulted 
in  behavior  of  the  viscous/inviscid  iteration  that  was  at 
variance  with  previous  experience.  The  boattail,  shown  in 
Figure  C-1,  consisted  of  a circular  arc  section  fairing 
smoothly  into  a conical  surface  of  15  degrees  which  then 
joined  the  cylindrical  sting  at  this  angle.  Shown  for  com- 
parison in  this  figure  is  the  NASA  Langley  = 1.768 

Circular  Arc  Boattail,  and  it  is  obvious  that  the  body  slopes 
of  the  Chow,  Bober,  Anderson  (CBA)  boattail  are,  for  a 
greater  distance,  steeper  than  the  NASA  body,  although  the 

closure  ratio  (d  /d  ) is  less, 
ex  m 

The  CBA  boattail  was  analyzed  at  four  Mach  numbers: 
0.56,  0.7,  0.8,  and  0.9.  The  data  used  for  comparison  were 
taken  from  the  experimental  work  of  Shrewsbury  (76).  As  can 
be  seen  in  Figures  C-2  and  C-3,  the  CVI  Program  predictions 
agree  well  with  the  data,  excepting  the  prediction  of  a higher 
pressure  near  the  boattail/sting  juncture,  which  has  been  en- 
countered previously  and  may  be  due  to  a small  separated 
region  of  an  inaccurate  value  of  6*.  At  Mach  number  0.8, 
however,  there  appears  to  be  an  instability  in  the  viscous/ 
inviscid  iteration,  in  that  two  solutions  are  produced. 

This  situation  is  best  shown  in 
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Figure  C-1.  Comparison  of  boattail  and  plume  simulator  shapes 
between  the  NASA  Langley  1.768  circular  arc  boattail  and 
the  CBA  boattail. 
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Figure  C-2.  Longitudinal  pressure  coefficient  distribution  for  the 
CBA  boattail,  M =0.56. 
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Figure  C-4,  in  which  the  repetitive  appearance  of  two  dis- 
similar solutions  appears  as  a large  difference  in 
from  one  iteration  to  the  next,  followed  by  a recovery  itera- 
tion and  an  apparently  converging  solution.  Relaxing  the 
boundary  layer  even  more  than  a factor  of  0.3  resulted  in  a 
longer  period  between  the  two  solutions,  but  little  else  was 
changed . 

The  two  solutions  are  shown  in  Figure  C-5.  It  is 
apparent  that  the  subsonic  solution,  the  more  stable  of  the 
two,  is  the  correct  solution,  while  the  supercritical  solu- 
tion, which  appears  once  at  intervals  depending  on  boundary 
layer  relaxation,  is  in  error,  particularly  in  the  expansion 
about  the  shoulder  of  the  boattail. 

At  a Mach  number  of  0.9,  the  solution  to  the  CVI 
Program  prediction  is  shown  in  Figure  C-6.  This  solution  con 
verged  quickly,  but  it  is  grossly  in  error.  It  does,  however 
protract  the  characteristic  of  the  erroneous,  or  "anomalous" 
solution  of  the  0.8  Mach  number  case,  i.e.  the  movement  of 

the  point  of  minimum  C downstream  one  grid  point. 

P 

At  the  time  of  this  writing,  the  reasons  for  such 
anomalous  behavior  are  still  in  question.  One  numerical 
experiment,  however,  in  which  the  inviscid  portion  of  the  CVI 
Program  was  restricted  to  central  differencing  is  worthy  of 
note.  In  the  0.8  Mach  number  case,  rotated  differencing  in 
the  South- Jameson  Transonic  Program  was  bypassed  after  itera- 
tion six.  The  solution  proceeded  to  converge  smoothly  and 


154 


Ln 


Figure  C*4.  Maximum  pressure  coefficient  change  at  each  iteration 
showing  presence  of  two  solutions,  = 0.80. 
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Figure  C-5.  Longitudinal  pressure  coefficient  distribution  for  the 
CBA  boattail  showing  presence  of  two  solutions,  = 0.80. 
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Figure  C-6.  Longitudinal  pressure  coefficient  distribution  for  the 
CBA  boattail  showing  erroneous  solution,  M - 0.90. 
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definitely,  and  the  resulting  C distribution  is  shown  in 

tr 

Figure  C-7.  It  agrees  well  with  the  data.  Similar  experi- 
mentation with  the  0.9  Mach  number  case  failed,  as  the  cen- 
tral differencing  was  not  capable  of  converging  under  super- 
sonic conditions  about  the  shoulder  of  the  boattail. 
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X,  inches 

Figure  C-7.  Longitudinal  pressure  coefficient  distribution  for  the 
CBA  boattail  using  central  differences  only,  = 0.80. 
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NOMENCLATURE 

Symbols  with  generally  recognized  definitions  are 
listed  as  such . Other  symbols  are  referred  to  the  applicable 
portion  of  the  text.  Additional  terms,  used  primarily  as  in- 
put parameters  to  particular  programs,  may  be  found  in  the 
applicable  table, 
a Local  speed  of  sound 

C Crocco  number  of  internal  flow  at  initial  con- 

a 

. . '^a 

ditions,  C_  = 

*^2 

P o 

Cj  Friction  coefficient,  Bartz , p. 30 

Cjj  Stanton  number,  Bartz,  p.  31 

C Pressure  coefficient 

P 

Cp^  Specific  heat  at  constant  pressure  for  internal 

flow 

d Diameter 

H H = 1 + Kn,  RAXBOD,  p.  19 

K Curvature  of  reference  coordinate  system, 

RAXBOD,  p.  19 
a Length 

M Local  Mach  number 

p Pressure 

Pr  Prandtl  number 

q Total  velocity , , RAXBOD,  p.  21 

r Radius  of  body,  Bartz,  p.  30 
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R Radius 

R Boundary  layer  relaxation  factor,  p.  44 

Re  Reynolds  number 

Re/Jl  Unit  Reynolds  number,  ft  ^ 

R Reynolds  number  based  on  momentum  thickness, 

0 

Bartz,  p.  31 

Reynolds  number  based  on  energy  thickness, 

'P 

Bartz,  p.  31 
T Temperature 

u,v  Streamwise  and  transverse  velocity  components, 

Chapman-Korst  mixing  theory,  p.  58 
U,V  streamwise  and  transverse  velocity  components 

U,V  Velocity  components  in  CiH  system,  RAXBOD , p.  19 

x,y  Intrinsic  coordinate  system,  Chapman-Korst 

mixing  theory,  p.  58 

X,Y  Streamwise  and  transverse  coordinates  in  the 

physical  coordinate  system 
2 Axial  coordinate,  Bartz,  p.  30 

6 Boundary  layer  thickness 

(S*  Boundary  layer  displacement  thickness 

e Eddy  viscosity,  Chapman-Korst  mixing  theory,  p.  58 

n Nondimens ional  similarity  variable,  Chapman-Korst 

mixing  theory,  p.  60 

Tip  Position  parameter,  Chapman-Korst  mixing  theory, 

p.  60 

0 Inclination  of  reference  coordinate  system, 

RAXBOD,  p.  19 
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6 Momentum  thickness,  Bartz,  p.  30 

A Stagnation  temperature  ratio,  Chapman -Korst  mixing 

theory,  p.  63 

5,n  -Orthogonal  curvilinear  coordinates,  RAXBOD,  p.  19, 

or  sheared  cylindrical  coordinates,  RAXBOD,  p.  20 
p Density 

0 Mixing  coefficient.  Chapman -Korst  mixing  theory, 

p.  58 

<J)  Disturbance  potential,  RAXBOD,  p.  19 

(|)  Energy  thickness,  Bartz,  p.  30 

^ Velocity  function , u/u_ , Chapman -Korst  mixing 

Q. 

theory,  p.  58 

4)  Nondimens ional  x-coordinate , x/6  , Chapman-Korst 

a 

mixing  theory,  p.  58 

SUBSCRIPTS 

a,b  Initial  internal  and  external  quantities,  Chapman- 

Korst  mixing  theory,  Figure  13,  p.  56 
aw  AdicODatic  wall 

bod  Body 

BT  Quantity  measured  from  beginning  of  boattail 

e Pertaining  to  effective  plume 

ex  Pertaining  to  exhaust 

j Pertaining  to  dividing  streamline  location, 

Chapman-Korst  mixing  theory.  Figure  13,  p.  56 
m Pertaining  to  intrinsic  coordinate  obtained  in 

Chapman-Korst  mixing  theory.  Figure  13,  p.  56 
m,max  Maximum 
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NN  Second  derivative  normal  to  velocity  vector, 

RAXBOD,  p.  21 

NOZ  Pertaining  to  nozzle 

o Stagnation  condition 

p Value  at  particular  x-locations  used  by  RAXBOD, 

p.  43 

R^,R^  Reference  limits  of  internal  and  external  flows 
rel  Pertaining  to  relaxed  boundary  layer  displacement 

thickness 

SS  Second  derivative  along  velocity  vector,  RAXBOD, 

p.  21 

w Wall  condition 

“ Free-stream  conditions 

2a, 2b  Initial  boundary  layer  velocity  profile  shapes 

I, II  One-stream  and  two-stream  values  of  a,  respectively,  „ 

Chapman-Korst  mixing  theory,  p.  79 

SUPERSCRIPTS 

’ First  derivative  in  longitudinal  direction 

" Second  derivative  in  longitudinal  direction 

n,n-l  Pertains  to  number  of  iteration  in  the  South- 


163 


AEDC-TH-77-106 


i2(n> 
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T 

T 
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2 A 

Ca  4. 


2 

b 
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